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SUMMARY 


New improved curve fits for the thermodynamic properties of 
equilibrium air have been developed. The curve fits are for 
p = p(e,p}» a = a(e,p) , T = T(e,p), s = s{e 5 p) 9 T = T(p,p)» 
h = h(p,p), p = p{p,s), e = e{p,s) and a = a(p,s). These curve 
fits can be readily incorporated into new or existing CFD codes 
if "real -gas effects" are desired. The curve fits were constructed 
using Grabau-type transition functions to model the thermodynami c 
surfaces in a piecewi se manner. The accuracies and continuity of 
these curve fits are substantially improved over those of previous 
curve fits appearing in NASA CR-2470. These improvements were due 
to the incorporation of a small number of additional terms in the 
approximating polynomial s and careful choices of the transition 
functions. The ranges of validity of the new curve fits are temper- 
atures up to 25,000 K and densities from 10“ 7 to 10 3 amagats (p/pq)- 


1 



SYMBOLS 


a 

a 0 
e 

e 0 
h 

P 

PO 

R 

s 

s O 

I 

T o 

Y 

P 

P 0 


speed of sound 

standard speed of sound = 331.36 m/s 

specific internal energy, m 2/s 2 

standard specific internal energy = 78,408.4 m ^/s 2 

specific enthalpy, m2/s2 

pressure, N/m 2 

standard pressure = 101,330.0 N/m^ 

gas constant = 287.06 m 2/s 2 - K 

specific entropy, m ^s 2 - K 

standard specific entropy = 6779.2 m2/s2 - K 

temperature, K 

standard temperature = 273. 15 K 
h/e 

density, kg/m 3 

standard density = 1. 292 kg/m^ 


Subscripts: 
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et al . (ref . 8) developed a 27-reaction equilibrium air program to 
generate the thermodynamic properties in tabular form as functions 
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not significantly increase the computat 



Using an approach of carefully choosing the Grabau-type transi- 
tion functions and using complete bicubic polynomials, the following 
curve fits were developed: 

P = p(e,p) 
a = a(e,p) 

T = T(e,p) 
s = s(e,p) 

T = T(p»p) 

h = h(p,p) 
p = p(p,s) 
e = e{p,s) 
a = a(p,s) 

These curve fits were based on the NASA RGAS data and have the same 
ranges of val idity, namely, temperatures up to 25,000 K and densities 
from ID"? to 10 ^ amagats (p/pq). 

BEHAVIOR OF AIR AT HIGH TEMPERATURE 

When a gas composed of polyatomic molecules is heated to high 
temperatures, its composition changes as a resul t of the chemical 
reactions which take place. Such a situation exists behind the 
shock wave which envelops a vehicle reentering the atmosphere of 
the earth. As a resul t of the change in chemical composition, the 
thermodynami c properties of the gas al so change. When the temperature 
of the gas is raised appreciably higher than the temperature at 
which dissociation reactions begin to occur, the electrons receive 
energy quanta because of the collisions between atoms. If the tem- 
perature , and hence the kinetic energy of the atoms is high enough 


7 



behavior provides a qualitative insight into the choice of the approx- 
imating functions. Figure 1 shows the function y plotted with respect 
to log 10 (p/p 0 ) - log 10 (p/p 0 ) at a density of 10" 7 amagats. Also 
shown are the various segments into which the curve may be divided, 
as indicated by A, AA, B, C and D. These segments are basically 
quadratic or 1 inear curves which are joined together by transition 



that electrons are removed from their orbits, ionization of the 
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curves. Two types of transition curves appear in figure 1 and these 
are illustrated in figures 2 and 3. Figure 2 shows a transition 
function which passes through a point of inflection and is referred 



Figure 3. A transition wi thout inflection 
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to as a transition with inflection. Figure 3 illustrates the second 
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Figure 4. Variation of y with log 10 (p/p 0 ) - log 10 (p/p 0 ) 


might be correlated with the simultaneous abrupt increases of the 
concentrations of ionized oxygen and nitrogen atoms and of ionized 
nitrogen molecules. 
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the terminology of Grabau (ref. 11) , are referred to as odd 
transitions. 

The second type of transition function is the kernel of the 
Bose-Einstein distribution function 


1 

1 - exp(kx) 


( 2 ) 
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Figure 6. Two odd transition functions 


whi ch provides transitions leading from one function to another 
wi thout a point of inflection and is obtained by merely changing 
the sign before the exponential term in the denominator of the Fermi - 
Dirac function. The transition function given by equation (2) is 
termed an even transition. Figure 7 illustrates two transitions 
of this kind between the x-axis and the 1 ine y = x, where, as before , 
the directions and rates of the transitions are governed by the sign 
and magnitude of the exponential constant k. It is important to 
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Figure 7 . Two even transition functions 


note that the expression for an even transition becomes an indeter- 
minate form when x is equal to the x-coordinate of the point of 
intersection of the two 1 ines bounding the transition. 

In the current study, each of the thermodynamic curves was 
approximated by means of quadratic or incomplete cubic segments 
connected by odd and even transitions as described above. Almost 
without exception, all the curves undergo odd transitions at low 
densities which gradual ly diminish as the density increases and 
change to even transitions. There are two ways of applying each 
of these transition functions. When the path of a curve appears 
to move from one straight 1 ine to another, there is present an offset 
which can be calculated in the direction of either of the variables. 
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From considerations of accuracy it appears to be better to view 
the transition in terms of the smaller offset. Both ways of viewing 
the offsets involve the choice of a base line. The use of the offset 
in the y-di recti on simpl ifies thi s choice since the x-axi s serves 
as a natural base line. 

Consider the problem of determining the equation of a curve 
con si sting of two 1 inear segments connected by an odd transition 
function (fig. 8) . The lower and upper 1 ine segments are given 
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However, they are found to be substantially alike in most instances. 


The determination of the constants of an even t ran si 
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is given by the ratio of the derivatives of the numerator and the 
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This approach of determining the constants of the Grabau- type 


transitions was extended in the present work to approximate transi- 
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transition function. The general technique of determining the values 
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where k is a positive constant ( typi cal ly, 20 < k < 25) chosen such 
that exp (k) and exp (-k) do not yield overflow or underflow conditions, 
respectively, on a computer. The constants a Q through a 3 can now 
be obtained in a straightforward manner from the system of four 
1 inear equations in four unknowns (eqs. ( 18)- ( 21 ) ) . 

The above method of obtaining the Grabau-type transition func- 
tions proved quite accurate in ensuring a negligible mi smatch in 
the dependent variable over the boundaries of adjoining sub-regions. 

It is a merit of this stepwi se method of constructing empirical equa- 
tions that any part can be removed for corrections without disturbing 
the surface approximation as a whole. 

EQUATIONS OF THE CURVE FITS 


The curve fits for the various thermodynami c properties were 
constructed using Grabau-type transition functions, as described 
previously. The general form of these curve fits can be written 
as 


z(x,y) = f^x.y) + 


f 2 (x,y) - f . (x,y) 

1 + exp{k Q + kjX + k^y + k^xy) 


( 22 ) 


where, in general 


f 1 (x,y) = Pj + P 2 x + p 3 y + p 4 xy + p g x + p g y + p y 


x y 


+ P 8 xy + p g x + p 1Q y 


(23) 
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Figure 9. Division of curve fit range by density 


dependent variables across these two division lines the following 
technique was adopted. If the choice of independent variables yielded 
a point within a specified band about either of these division lines, 
the dependent variable was linearly interpolated between the values 
obtained at the end points of the band. 

P = p(e,p) 

For the correlation of p = p(e,p ) , the ratio y = h/e was curve- 
fitted as a function of e and p so that p can be calculated from 
p = pe(y - 1) (25) 


25 




The genera] form of the equation used for y was 

Y = a x + a £ Y + a 3 Z + a 4 YZ + a 5 Y 2 + a g Z 2 + a ? Y 2 Z + a g YZ 2 + a g Y 3 

+ a 1 Q Z 3 + (a n + a 12 Y + a 13 Z + a^YZ + a 15 Y 2 + a lg Z + a iy Y Z 
+ a ^YZ t a^gY + a 20 ^ )/L 1- oxp{a 2 i "*• ®22^ ^23^ ®24^'^ 

(26) 

where Y = log J p / p J and Z = log fe/RT ). The units for p are kg/m 3 

10 0 10 v 

and the units for e are m 2 /s 2 . It should be noted that some of 

the terms appearing in the above equation are not used over the 

complete range of e and p . 

a = a(e,p) 

An exact expression for the speed of sound in terms of y was 
deri ved by Barnwel 1 (ref. 13) and may be wri tten as 


a 


e (y - 1) 


Y + 


8Y 


3log e e 


3? 


Slog p 


1 1/2 


(27) 


Since complete bicubic polynomials were used for f j(Y,Z) and 
f 2(Y,Z) - f j(Y,Z) in equation (26) for y, equation (27) was used 
directly for the correlation a = a(e,p) without further corrections, 
uni ike references 10 and 15. 
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T = T(e,p) 

In the calculation of T = T(e,p), the pressure is first deter- 
mined using equation (25) , and then the temperature is calculated 
using the equation 

1og 10 (T/T 0 ) = b 2 + b 2 Y + b 3 Z + b 4 YZ + b 5 Y 2 + b 6 Z 2 + b ? Y 2 Z + bgYZ 2 

+ V 3 + b 10 z3 + (b ll + b 12 Y + b 13 Z + b 14 YZ + b 15 y2 

+ ^15^ + b i7^ ^ + ^l8 yZ + b l9 Y + ^20 Z + ex P(^21 

+ b 22 Y + b 23 Z + b 24 YZ)] 

(28) 


where Y = log 10 (p/p 0 ), X = log 10 (p/p 0 ), and Z = X - Y. The units 
for p are newtons/m , and the units for T are kelvin. The coeffi- 
cients bj to b 24 were determined in such a way as to compensate 
for the errors incurred in the initial calculation of pressure using 
equation (25). 

h = h(p >P ) 

For the correlation of h = h(p,p), the ratio y = h/e was curve fitted 
as a function of p and p so that h can be calculated from 

h * (p/p) & (29) 

The general form of the equation used for y was 
y = Cj + c 2 Y + c 3 Z + c 4 YZ + c 5 Y 2 + Cg Z 2 + c y Y 2 Z + Cg YZ 2 + c g Y 3 

+ c 10 z3 + (c ll + C 12 V + C 13 Z + C 14 YZ + c 15 y2 + c 16 z2 + c 17 y2z 
+ c 18 YZ 2 * c 19 Y 3 * c 20 Z 3 )/[1 + exp(c 21 + c 22 Y + c 23 Z + c 24 YZ)] 

( 30 ) 
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where Y = log 10 ( P / P() ), X = log 10 (p/p Q ), and Z = X - Y. For the 
correlations p = p(e,p) and h = h(p,p)> where y was the variable 
curve-fitted, an even transition function was used to model the transi- 
tion between the perfect gas equation and the remainder of the curve 
fit in the lowest density region (- 7.0 < log 10 (p/p Q ) < - 4.50). 

This yielded a more accurate fit than an ordinary bi cubi c curve 
without any transitions, though the latter seems intuitively more 
appropriate. 

T = T(p,p) 

The general form of the equation used for the correlation 
T = T(p,p) was 

log 10 (T/T 0 ) = dj + d 2 Y + d 3 Z + d 4 YZ + d,-Y 2 + d 6 Z 2 + d ? Y 2 Z + dgYZ 2 
♦ d,Y 3 ♦ d w Z 3 ♦ <d u ♦ d 12 Y + d 13 Z ♦ d^YZ ♦ d 16 Y 2 
+ d 16 Z 2 + d 17 Y 2 Z + d 18 YZ 2 + d 1? Y 3 + d 2Q Z 3 )/[l + exp(d 2J 



( 32 ) 
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where Y = log^Cp/p^ and Z = log 10 (e/RT 0 ) . The units for s are 
m^/s^ - K. As is evident from equation (32), Grabau transition 
functions were not necessary for thi s curve fit. 

p = p(p,s) 

Uni ike the preceding curve fits in which density is one of 
the independent variables, the domain of the curve fit p = p (p,s) , 
as well as the curve fits e = e(p,s) , and a = a(p,s) , cannot be 
divided into subdomains on the basis of density. For reasons of 
accuracy, it was necessary to subdivide the domain in terms of s 
as seen in figure 10. 



L0G 10 (p/p 0 ) 


Figure 10. Divi si on of curve fit range by entropy 
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The general form of the equation used for the correlation 
p = p (p,s) was 

1 °gi0 ( p /p 0 ) = f i + f 2 y + f 3 Z + f 4 YZ + f 5 y2 + f 6 Z 2 + f 7 y2z + yz 2 

+ V 3 + f 10 z3 + (f ll + f 12 v + f 13 z + f 14 YZ + f 15 y2 
+ V 2 + f 17 Y2z + f 18 Vz2 + f 19 y3 + f 20 z3 > / f 1 + ex P< f 21 
+ f 22 Y + f 23 Z ^ 

(33) 

where Y = logy s/R) , X = logy p/p g ) , and Z = X - Y. The units 
for s are m z /sec z - K. The denominator of the Grabau transition 
function in equation ( 33) is 1 inear in Y and Z and not quadratic. 

e = e(p,s) 

For the correlation of e = e(p,s) , the general form of the curve-fit 
equation was 

log 10 (e/RT 0 ) = gj + g 2 Y + g 3 Z + g 4 YZ + g 5 Y 2 + g g Z 2 + g y Y 2 Z + g g YZ 2 

+ g 9 V3 + 9 10 Z3 * (g ll + 9 12 Y + g 13 Z + 9 14 YZ + g 15 y2 

+ g 16 Z 2 + g 17 Y 2 Z + g 18 YZ 2 + g ig Y 3 + g 2Q Z 3 )/[l + exp(g 21 

+ g 22 Y + 9 23 z)] 

(34) 
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where Y = log^s/R), X = log^Cp/Pg), and Z = X - Y. 

a = a(p,s) 

For the correlation of a = a(p,s), the general form of the equation 
used was 

log 10 (a/a 0 ) = hj + h 2 V * hjZ + h 4 YZ + h/ + h 6 Z 2 + h ? Y 2 Z + h 8 YZ 2 

+ V 3 + h 10 Z3 + (h ll * h 12 Y + h 13 Z + h 14 YZ + h 15 y2 
+ hjgZ 2 + h 1? Y 2 Z + hjgYZ 2 + h ig Y 3 + h 2 „Z 3 )/[l + exp(h n 

+ h 22 Y + h 23 Z)] 

(35) 


where Y = log^s/R), X = log 1() (p/p 0 ), and Z = X - Y. The units 
of a are m/sec. 

RESULTS AND CONCLUSIONS 

New simplified curve fits for the thermodynami c properties 
of equilibrium air were constructed using the procedures described 
in the preceding chapters. Comparisons of the curve fits p = p(e,p) , 
a = a(e,p), T = T(e,p), s = s(e,p), T = T(p,p), h = h(p,p),p =p(p,s), 
e = e(p,s) and a = a(p,s) with the original NASA RGAS program are 
shown in figures 11 to 19. The following procedure was empl oyed 
in making the comparisons for the first four curve fits. First, 
p and p data were suppl ied as input to the NASA RGAS program and 
e was computed. Then, this e and the original p were inputed into 
the TGAS1 subroutine to obtain p, a and T and into the TGAS2 sub- 
routine to obtain s. As a result of this procedure, log 1Q (p/p 0 ) 
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NASA RGAS 



igure 12. Comparison of 
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at log in (p/p 



TABLE 2. ACCURACY OF TGAS1 FOR a = a(e,p) 
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and T = 1.25 



at log 10 (p/pg) = - 0.625; log^e/RTg) = 3.255 
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at log in ( p /p n ) = - 0.625; log in (p/p n ) - 1og in (p/p n ) = 2.44 



TABLE 6. ACCURACY OF TGAS4 FOR h = h(p, P ) 
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at 1og.Jp/p 



TABLE 7. ACCURACY OF TGAS5 FOR 
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TABLE 9. ACCURACY OF TGAS7 FOR a = a(p,s) 
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Figure 20. Example curve fit for p = p(e,p ) 


(ref. 15) for the same curve fits is 3.4 times faster than the 
NASA RGAS subroutine. The new TGAS3 subroutine for T = T(p,p) is 
2.7 times faster than the NASA RGAS subroutine as compared to the 



TABLE 10. COMPARISON OF VARIABLES AT JUNCTURE POINTS FOR THE CURVE 


FIT p = p(e,p) 


„ . POINT A POINT B POINT C POINT D POINT E 

Density 


ratio 

p/p 0 

lower 

upper 

lower 

upper 

lower 

upper 

lower 

upper 

lower 

upper 

o 

i 

"vi 

1.79xl0' 2 

1.81xl0" 2 

7.32xl0' 2 

7. 38xl0' 2 

1. 90x10" 1 

1.90X10' 1 

8.72X10' 1 

8.72X10" 1 

2.62x10° 

2.63x10° 

in " 6 

1.79X10' 1 

1. 80x10' 1 

7.81x10'“ 

7. 78x10" 1 

2.01x10° 

2.03x10° 

9.53x10° 

9.63x10° 

2.86X10 1 

2.89X10 1 

io " 5 

1.79x10° 

1.80x10° 

8.17x10° 

8.19x10° 

2.15X1C 1 

2.18X10 1 

1.05xl0 2 

1.06xlQ 2 

3. 15xl0 2 

3. 16xl0 2 

K)" 4 

1.80X10 1 

1.81X10 1 

8.67X10 1 

8.70X10 1 

2.40xl0 2 

2.43xl0 2 

9.78xl0 2 

9.79xl0 2 

1.80xl0 2 

1.81xl0 2 

io " 3 

1.80xl0 2 

1.81xl0 2 

9.12xl0 2 

9. 13xl0 2 

2.61xl0 3 

2.63X10 3 

1.09xl0 4 

1.09xl0 4 



io " 2 

1.80xl0 3 

1.81xl0 3 

9. 51xl0 3 

9. 51xl0 3 

2.83xl0 4 

2.84xl0 4 

1.23xlO S 

1.23xlO S 



10" 1 

1.80xl0 4 

1.81xl0 4 

9.80xl0 4 

9.81xl0 4 

3.08xl0 5 

3.08xl0 5 

1.39xlO S 

1.39xlO S 



10° 

1.80xl0 5 

1.81xl0 5 

1. 36xl0 6 

1. 36xl0 6 

4.11xl0 6 

4.11x10 s 





IQ 1 

1.80xl0 6 

1.81xl0 6 

1.41xl0 7 

1.41xl0 7 

4.50xl0 7 

4.54xl0 7 





10 2 

1.80xi0 7 

1.81xl0 7 

1.45x10® 

1.43xl0 8 

4.91xl0 8 

5.00x10® 





10 3 

1.80xl0 8 

1.83xl0 8 

1.48xl0 9 

1.45xl0 9 

5.42xl0 9 

5.53xl0 9 
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TABLE 11. COMPARISON OF VARIABLES AT JUNCTURE POINTS FOR THE CURVE 
FIT a = a (e ,p ) 


Density 

ratio 

p/p 0 

POINT A 

POINT B 

POINT C 

POINT D 

POINT E 

lower 

upper 

lower 

upper 

lower 

upper 

lower 

upper 

lower 

upper 

10- 7 

440 

441 

769 

790 

1250 

1260 

2718 

2733 

4731 

4715 

10- 6 

440 

439 

808 

814 

1291 

1307 

2857 

2871 

4983 

5016 

10- 5 

440 

438 

831 

841 

1343 

1359 

3021 

3029 

5259 

5287 

10- 4 

441 

440 

869 

874 

1429 

1441 

2923 

2925 



10' 3 

441 

440 

902 

904 

1498 

1506 

3115 

3116 



10" 2 

441 

440 

932 

932 

1573 

1578 

3337 

3341 



10' 1 

441 

441 

957 

957 

1655 

1656 

3596 

3602 



10° 

442 

441 

1120 

1118 

1924 

1924 





10 1 

442 

440 

1149 

1145 

2027 

2039 





10 2 

442 

440 

1171 

1164 

2141 

2166 





10 3 

442 

441 

1188 

1179 

2287 

2312 
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TABLE 12. COMPARISON OF VARIABLES AT JUNCTURE POINTS FOR THE CURVE 
FIT T = T(e,p ) 


POINT a point b point c point d point e 

Density 

ratio ___ 

p/p 0 lower upper lower upper lower upper lower upper lower upper 


1G~ 7 

486 

481 

2112 

2091 

4033 

4034 

7868 

7869 

10- 6 

486 

482 

2181 

2168 

4283 

4284 

8471 

8479 

10- 5 

486 

484 

2243 

2243 

4548 

4548 

9145 

9146 

IQ" 4 

486 

481 

2312 

2312 

4837 

4818 

10,364 

10,319 

Kf 3 

486 

481 

2347 

2366 

5090 

5088 

11,190 

11,177 

10' 2 

486 

481 

2376 

2404 

5307 

5326 

11,958 

12,006 

10" 1 

486 

481 

2400 

2417 

5508 

5517 

12,702 

12,738 

10° 

486 

482 

2408 

2414 

6242 

6265 



10 1 

486 

482 

2413 

2416 

6585 

6595 



10 2 

486 

482 

2416 

2416 

6955 

6960 



10 3 

486 

483 

2418 

2419 

7317 

7328 
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TABLE 13. COMPARISON OF VARIABLES AT JUNCTURE POINTS FOR THE CURVE 
FIT T = T(p,p ) 


Den si ty 
ratio 
p/p 0 

POINT 

A 

POINT 

B 

POINT 

C 

POINT 

D 

POINT E 

lower 

upper 

lower 

upper 

lower 

upper 

lower 

upper 

lower upper 

10- 7 

486 

482 

2089 

2089 

4025 

4033 

7864 

7838 


10- 6 

486 

482 

2165 

2165 

4281 

4281 

8470 

8481 


-5 

10 3 

486 

484 

2242 

2242 

4549 

4554 

9146 

9146 


10" 4 

486 

482 

2310 

2310 

5064 

5042 

10,796 

10,746 


10" 3 

486 

481 

2363 

2363 

5386 

5376 

11,793 

11,682 


10” 2 

486 

481 

2404 

2404 

5690 

5701 

12,742 

12,679 


10" 1 

486 

482 

2402 

2402 

5968 

5998 

13,671 

13,687 


10° 

486 

482 

2700 

2700 

6248 

6267 




10 1 

486 

482 

2706 

2710 

6585 

6598 




IQ 2 

486 

483 

2711 

2712 

6950 

6959 




10 3 

486 

483 

2713 

2713 

7309 

7319 
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TABLE 14. COMPARISON OF VARIABLES AT JUNCTURE POINTS FOR THE CURVE 
FIT h = h(p,p ) 


POINT A POINT B POINT C POINT D POINT E 

Density 

ratio — — ~ 

p/Pq lower upper lower upper sower upper lower upper lower upper 


10" 7 

0.346xl0 6 

0. 346xl0 6 

0.282xl0 7 

0. 285xl0 7 

0.160xl0 8 

0.159X10 8 

0.997xl0 8 

0.997xl0 8 

!0' 6 

0. 346xl0 6 

0.346xl0 6 

0.253xl0 ? 

Q.254xlQ 7 

0. 138xl0 8 

O 

0.1 

CO 

X 

o 

00 

0.890xl0 8 

0.890xl0 8 

10- 5 

0. 346xl0 6 

0.346xl0 6 

0.233xl0 7 

0.235xl0 7 

0.120xl0 8 

0.122xl0 8 

0.793xl0 8 

0.792xl0 8 

10" 4 

0.346xl0 6 

0. 346xlQ 6 

0.345xl0 7 

0. 345xl0 7 

0.247xl0 8 

0.247xl0 8 

0.812xl0 8 

0.813xl0 8 

10- 3 

0.346xl0 6 

0.346xl0 6 

0. 3 1 4x 10 7 

0. 315xl0 7 

0. 214xl0 8 

00 

o 

X 

C\j 

o 

0.720xl0 8 

0.721xl0 8 

10" 2 

0. 346xl0 6 

0.346xl0 6 

0. 296xl0 7 

0.296xl0 7 

00 

o 

X 

wo 

CO 

o 

0. 186xl0 8 

0.646xl0 S 

0.646xl0 8 

10“ 1 

0. 346xi0 6 

0.346xl0 6 

0.288xl0 7 

0.288xl0 7 

0.164X10 8 

0. 164xl0 8 

0. 590xl0 8 

0. 591xl0 8 

10° 

0.345xl0 6 

0.345xl0 6 

0. 386xl0 7 

0. 387 xlO 7 

o 

r\3 

o 

X 

o 

00 

0. 202xl0 8 



10 1 

0. 345xl0 6 

0. 345xl0 6 

0.377xl0 7 

0. 38QxlQ 7 

00 

o 

X 

o 

00 

o 

0.181xl0 8 



10 2 

0.345xl0 6 

0.345xl0 6 

0.374xl0 7 

0.376xl0 7 

0.166xlQ 8 

0.166xl0 8 



10 3 

0.345xl0 6 

0. 345xlQ 6 

0. 374xl0 7 

0.374xl0 7 

0. 156xl0 8 

0.156xl0 8 




57 



TABLE 15. COMPARISON OF COMPUTER TIMES 


Curve Fit 


p = p(e,p) 
a = a (e , p) 
T = T(e, p) 

s = s(e, p ) 

T = T(p, p) 

h = h(p, p) 

p = p(p,s) 
e = e(p,s) 
a = a(p,s) 


Number of 
Data Points 


Computer Time (seconds) 


Old 


New 


Subroutine Subroutine 


NASA RGAS 


10,661 

0.54 

0.77 

1.86 

10,661 

— 

0.20 

2.03 

9,921 

0.25 

0.31 

0.84 

9,921 

0.19 

0.26 

0.84 

3,038 

— 

0.10 

1.07 

3,038 

— 

0.11 

1,06 

3,038 


0.11 

1.06 
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previous TGAS subroutine which is 3.4 times faster. The new TGAS4 
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APPENDIX 


All the curve fits developed in this study have been incorpor- 
ated into a single master program called TGAS. This master program 
permits the user to select the desired curve fit(s) from a "menu" 
of possibilities. The calling statement for this subroutine is 
CALL TGAS (P,RHO,E,H,T,A,S, NTGAS) 

with 

2 

P = pressure, N/m 

3 

RHO = density, kg/m 
E = specific internal energy, m^/s^ 

H = specific enthalpy, m ^/s^ 

T = temperature, K 

A = speed of sound, m/s 

S = specific entropy, m^/s^ - K 

NTGAS = an integer flag to be set by the user for selection 
of the appropriate curve fit(s) , 

NTGAS =0: p = p(e,p) 

NTGAS = lr p = p(e,p) , a = a (e ,p ) 

NTGAS =2: p = p(e,p), T = T(e,p) 

NTGAS =3: p = p(e,p), a = a(e,p), T = T(e,p) 

NTGAS =4: s = s(e,p) 

NTGAS =5: T = T(p,p) 

NTGAS =6: h = h(p,p) 

NTGAS =7: P = p(p,s) 

NTGAS =8: e = e(p,s) 

NTGAS =9: a = a(p,s) 
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The curve fits for p = p(e, p) , a = a(e, p) and T = T(e, p) have 
been placed in a single subroutine TGAS1. Subroutine TGAS2 computes 
s = s(e, p) , subroutine TGAS3 computes T = T(p, p) , and subroutine 
TGAS4 computes h = h(p,p). The curve fits p = p(p,s), e = e(p,s) 
and a = a(p,s) have been placed in subroutines TGAS5, TGAS6 and 
TGAS7 , respectively. The subroutines TGAS1 to TGAS7 can be used 
in a "stand-alone" manner if so desired, independent of the master 
program. A FORTRAN 1 i sting of these subroutines follows. 
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SUBROUTINE TGAS(P,R,E,H,T,A,S,NTGAS) 
PARAMETERS : 

P = PRESSURE, IN NEWTONS/M**2 
R = DENSITY, IN KG/M**3 
E = SPECIFIC INTERNAL ENERGY, IN (M/SEC)**2 
H = SPECIFIC ENTHALPY, IN (M/SEC)**2 
T = TEMPERATURE, IN KELVIN 
A = SPEED OF SOUND, IN M/SEC 
S = SPECIFIC ENTROPY, IN (M/SEC)**2/KELVIN 
NTGAS = INTEGER FLAG TO BE SET AS FOLLOWS , 


NTGAS = 0 : 

INPUT : R,E ; OUTPUT : 

NTGAS = 1 : 

INPUT : R,E ; OUTPUT 

NTGAS = 2 : 

INPUT : R,E ; OUTPUT 

NTGAS = 3 : 

INPUT : R,E ; OUTPUT 

NTGAS = 4 : 

INPUT : R,E ; OUTPUT 

NTGAS = 5 : 

INPUT : P,R ; OUTPUT 

NTGAS = 6 : 

INPUT : P,R ; OUTPUT 

NTGAS = 7 : 

INPUT : P,S ; OUTPUT 

NTGAS = 8 : 

INPUT : P,S ; OUTPUT 

NTGAS = 9 : 

INPUT : P,S ; OUTPUT 

IF ( NTGAS. EQ.O) CALL 
IF ( NTGAS. EQ.l) CALL 
IF (NTGAS.EQ.2) CALL 
IF (NTGAS. EQ. 3) CALL 
IF (NTGAS.EQ.4) CALL 
IF (NTGAS.EQ.5) CALL 
IF (NTGAS. EQ. 6) CALL 


P 


P.A 


P,T 


: S 


: T 


: H 


: R 


: E 


: A 

TGAS1(E,R,P,A,T,0) 

TGAS1(E,R,P,A,T,1) 

TGAS1(E,R,P,A,T,2) 

TGAS1(E,R,P,A,T,3) 

TGAS2(E,R,S) 

TGAS3(P,R,T) 

TGAS4(P,R,H) 
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IF (NTGAS.EQ.7) CALL TGAS5(P,S,R) 
IF (NTGAS.EQ.8) CALL TGA$6(P,S,E) 
IF (NTGAS.EQ.9) CALL TGAS7(P,S,A) 
RETURN 


KO 

kO 


Q 

Z 

UJ 



O Cl O O o o o oooooooooo 


SUBROUTINE TGAS1 (E,R,P,A J.MFLAG) 

INPUTS FOR SUBROUTINE: 

E= INTERNAL ENERGY, IN (M/SEC)**2 ; 

R= DENSITY, IN KG/M**3 . 

OUTPUTS: 

P= PRESSURE, IN NEWT0NS/M**2 ; 

A= SPEED OF SOUND, IN M/SEC ; 

T= TEMPERATURE, IN KELVIN . 

IF MFLAG=0, RETURN P ; 

IF MFLAG=1 , RETURN P AND A ; 

IF MFLAG=2 , RETURN P AND T ; 

IF MFLAG=3, RETURN P, A AND T . 

DATA EO.RO, PO, TO, GASCON/78408 . 4E00 , 1 . 292E00 , 1 . 0133E05 , 
*273 . 15E00 , 287 . 06E00/ 

RRATIO=R/RO 
ERATIO=E/EO 
Y=AL0G10( RRATIO) 

Z=AL0G10( ERATIO) 

LFLAG=0 
KFLAG=0 

IF (MFLAG.GT. 1) LFLAG=1 
IF ((MFLAG.EQ. 1) .OR.(MFLAG.EQ.3)) KFLAG=1 
I F(ABS( Y+4 . 5E00) . LT . 2 . 5E-02) GO TO 20 
I F( ABS( Y+0 . 5E00) . LT . 5 . 0E-03) GO TO 50 
IFLAG=-1 
GO TO 90 
10 IF(LFLAG.EQ.l) GO TO 300 

RETURN 

20 IFLAG=0 

RSAVE=R 
YM=Y 

Y=-4 . 5E00+2 . 5E-02 
YHIGH=Y 
R=( 10 . **Y)*R0 
JFLAG=-1 
GO TO 90 
30 PHIGH=P 

AHIGH=A 

Y=-4 . 5E00-2 . 5E-02 
YLOW=Y 

R=(10.**Y)*R0 
JFLAG=0 
GO TO 90 
40 PLOW=P 
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ALOW=A 
GO TO 80 
50 IFLAG=1 

RSAVE=R 
YM=Y 

Y=-0 . 5E00+0 . 5E-02 
YHIGH=Y 
R=( 10 . **Y)*RO 
JFLAG=-1 
GO TO 90 
60 PHIGH=P 

AHIGH=A 

Y=-0 . 5E00-0 . 5E-02 
YLOW=Y 

R=( 10 . **Y)*R0 
JFLAG=0 
GO TO 90 
70 PLOW=P 

ALOW=A 

80 P=PLOW+(PHIGH-PLOW)/(YHIGH-YLOW)*(YM-YLOW) 

A=A LOW+ ( AH I GH- A LOW ) / ( YH I GH- Y LOW ) * ( YM- Y LOW ) 

R=RSAVE 

IF(LFLAG.EQ.l) GO TO 300 
RETURN 

90 CONTINUE 

IF(Y.GT.-O.SEOO) GO TO 200 
IF(Y.GT.-4.5E00) GO TO 150 
IF(Z.GT.0.65E00) GO TO 100 
GAMM=1 . 3965E00 
GO TO 250 

100 IF(Z.GT. 1.5E00) GO TO 110 

GAS1=1 . 52792E00-1 . 26953E~02*Y 

GAS2=(-6.13514E-01-5.08262E-02*Y)*Z 

GAS3=(-5 . 49384E-03+4 . 75120E-05*Z-3 . 18468E-04*Y)*Y*Y 

GAS4=(6 . 31835E-01+3 . 34012E-02*Y-2 . 19921E-01*Z)*Z*Z 

GAS5=-4 . 96286E01— 1 . 17932E+01*Y 

GAS6=(6 . 91028E01+4 . 40405E+01*Y)*Z 

GAS7=( 5 . 09249E00-1 . 40326E00*Z+2 . 08988E-01*Y)*Y*Y 

GAS8=( 1 . 37308E01-1 . 78726E01*Y-1 .86943E01*Z)*Z*Z 

GAS9=EXP( 24 . 60452E00-2 . E00*Y-2 . 093022E01*Z) 

DEN0=1 . -GAS9 

GAMM=GAS 1+GAS2+GAS 3+GAS4+( GAS 5+GAS6+GAS7+GAS8 ) /DENO 
IF(KFLAG.EQ.O) GO TO 260 
GAS1R=-1 . 26953E-02 
GAS2R=-5.08262E-02*Z 

GAS'3R=( -1 . 098768E-02-9 . 50240E-05*Z-9 . 554040E-04*Y)*Y 

GAS4R=3.34012E-02*Z*Z 

GAS5R=-1 . 17932E01 

GAS6R=4. 40405E01*Z 

GAS7R=( 1 .018498E01-2 . 80652E00*Z+6 . 269641E-01*Y)*Y 
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GAS8R=-1.78726E01*Z*Z 

GAS9R=-2.0E00 

GAS2E=GAS2/Z 

GAS3E=4.75120E-05*Y*Y 

GAS4E=( 1 . 26367E00+6 . 68024E-02*Y-6 . 59763E-01*Z)*Z 

GAS6E=GAS6/Z 

GAS7E=-1 . 40326E00*Y*Y 

GAS8E=(2. 74616E01-3 . 57452E01*Y-5 . 60829E01*Z)*Z 
GAS9E=-2 . 093022E01 

GAMMR=GAS1R+GAS2R+GAS3R+GAS4R+(GAS5R+GAS6R+GAS7R+GAS8R)/DEN0 

*+(GAS5+GAS6+GAS7+GAS8)*GAS9R*GAS9/(DEN0**2) 

GAMME=GAS2E+GAS3E+GAS4E+(GAS6E+GAS7E+GAS8E)/DEN0 

*+(GAS5+GAS6+GAS7+GAS8)*GAS9E*GAS9/(DEN0**2) 

GO TO 260 

110 IF(Z.GT.2.2E00) GO TO 120 

GAS1=-1 . 70333E01-5 . 08545E-01*Y 
GAS2=(2.46299E01+4.45617E-01*Y)*Z 
GAS3=(-8.95298E-03+2.29618E-03*Z-2.89186E-04*Y)*Y*Y 
GAS4=(-1 . 10204E01-9 .89727E-02*Y+1 .62903E00*Z)*Z*Z 
GAS 5=1 . 86797E01+5 . 19662E-01*Y 
GAS6=(-2 .41338E01-4. 34837E~01*Y)*Z 
GAS7=( 9 . 16089E-03-1 . 52082E-03*Z+3 . 46482E-04*Y)*Y*Y 
GAS8=( 1 . 02035E01+9 . 70762E-02*Y-1 . 39460E00*Z)*Z*Z 
GAS9=(-1 . 42762E02-1 . 647088E00*Y+7 . 660312E01*Z 
*+8.259346E-01*Y*Z) 

IF(KFLAG.EQ.O) GO TO 240 
GASlR=-5 . 08545E-01 
GAS2R=4.45617E-01*Z 

GAS3R=(-1 . 790596E-02+4. 59236E-03*Z-8. 67558E-04*Y)*Y 

GAS4R=-9.89727E-02*Z*Z 

GAS5R=5 . 19662E— 01 

GAS6R=-4.34837E-01*Z 

GAS7R=( 1.832178E-02-3 . 04164E-03*Z+1 . 039446E-03*Y)*Y 

GAS8R=9 . 70762E-02*Z*Z 

GAS9R=-1 . 647088E00+8 . 259346E-01*Z 

GAS2E=GAS2/Z 

GAS3E=2 . 29618E~03*Y*Y 

GAS4E=( -2 . 20408E01-1 . 979454E-01*Y+4 . 88709E00*Z)*Z 

GAS6E=GAS6/Z 

GAS7E=-1.52082E-03*Y*Y 

GAS8E=( 2 . 0407E01+1 . 941524E-01*Y-4 . 1838E00*Z)*Z 
GAS9E=7 . 660312E01+8 . 259346E-01*Y 
GO TO 240 

120 IF (Z.GT.3.05E00) GO TO 130 
GAS 1=2 . 24374E00+1 . 03073E-01*Y 
GAS2=(-5 . 32238E-01-5 . 59852E-02*Y)*Z 
GAS3=( 3 . 56484E-03-1 . 01359E-Q4*Z+1 . 59127E-04*Y)*Y*Y 
GAS4=( -4 . 80156E-02+1 . 06794E-02*Y+3 . 66035E-02*Z)*Z*Z 
GAS5=-5 . 70378E00-3 . 10056E-01*Y 
GAS6=(5 . 01094E00+1 .80411E-01*Y)*Z 
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GAS7=( —9 . 49361E-03+1 . 94839E-03*Z-2 . 24908E-04*Y)*Y*Y 
GAS8=(-1 .40331E00-2.79718E-02*Y+1 .20278E-01*Z)*Z*Z 
GAS9=( 1 . 139755E02-4 . 985467E00*Y-4 . 223833E01*Z 
*+2 . 009706E00*Y*Z) 

IF(KFLAG.EQ.O) GO TO 240 
GAS1R=1 .03073E-01 
GAS2R=-5 . 59852E-02*Z 

GAS3R=(7 . 12968E-03-2 . 0218E-04*Z+4. 77381E-04*Y)*Y 
GAS4R=1 . 06794E-Q2*Z*Z 
GAS5R=-3 . 10056E-01 
GAS6R=1.80411E-01*Z 

GAS7R=(-1 . 898722E-02+3 . 89678E-03*Z-6 . 74724E-04*Y)*Y 

GAS8R=-2.79718E-02*Z*Z 

GAS9R=-4 . 985467E00+2 . 009706E00*Z 

GAS2E=GAS2/Z 

GAS3E=-1.01359E-04*Y*Y 

GAS4E=(-9 .60312E-02+2 . 13588E-02*Y+1 .098105E-01*Z)*Z 

GAS6E=GAS6/Z 

GAS7E=1.94839E-03*Y*Y 

GAS8E=( -2 . 80662E00-5 . 59436E-02*Y+3 . 60834E-01*Z)*Z 
GAS9E=-4 . 223833E01+2 . G09706E00*Y 
GO TO 240 

130 IF(Z.GT.3.4E00) GO TO 140 

GASl=-0 . 20807E02+0 . 4Q197E00*Y 
GAS2=(0. 22591E02-0 .25660E00*Y)*Z 
GAS3=(-0 . 95833E-03+0 . 23966E-02*Z+0 . 33671E-03*Y)*Y*Y 
GAS4=(-0 . 77174E01+0 . 4606E-01*Y+0 . 878E00*Z)*Z*Z 
GAS5=-0 . 21737E03-0 . 46927E01*Y 
GAS6=(0. 18101E03+0 . 26621 E01*Y)*Z 
GAS7=(-0 . 34759E-01+0 . 64681E-02*Z-0 . 70391E-03*Y)*Y*Y 
GAS8=(-0 . 50019E02-0. 38381E00*Y+0 . 45795E01*Z)*Z*Z 
GAS9=(0 . 4544373E03+0 . 1250133E02*Y-0 . 1376001E03*Z 
*-0 . 3641774E01*Y*Z) 

IF(KFLAG.EQ.O) GO TO 240 
GAS1R=0 . 40197E00 
GAS2R=-0 . 25660E00*Z 

GAS3R=(-1.91666E-03+4.7932E-03*Z+1.01013E-03*Y)*Y 

GAS4R=0.4606E-01*Z*Z 

GA$5R=-0 . 46927E01 

GAS6R=0 . 26621 EO 1*Z 

GAS7R=(-6 . 9518E-02+1 . 29362E-02*Z-2 . 11173E-03*Y)*Y 

GAS8R=-0. 38381E00*Z*Z 

GAS9R=0 . 1250133E02-0 . 3641774E01*Z 

GAS2E=GAS2/Z 

GAS3E=0.23966E-02*Y*Y 

GAS4E=(-1 . 54348E01+9 .212E-02*Y+2 . 634E00*Z)*Z 

GAS6E=GAS6/Z 

GAS7E=0.64681E-02*Y*Y 

GAS8E=(-1 . 00038E02-7 . 6762E-01*Y+1 . 37385E01*Z)*Z 
GAS9E=-0 . 1376001 E03— 0 . 3641774E01*Y 
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GO 10 240 

140 IF(Z.GT.3.69E00) WRITE(6 , 1000) R,E 

GAS l=-5.22951E01-4. 0001 1E-01*Y 
GAS2=(4. 56439E01+2 . 24484E-01*Y)*Z 
GAS3=(-3 . 73775E-03+2 . 43161E-03*Z+2 . 24755E-04*Y)*Y*Y 
GAS4=(-1 . 29756E01— 2 . 79517E-02*Y+1 . 22998E00*Z)*Z*Z 
GAMM=GAS1+GAS2+GAS3+GAS4 
IF(KFLAG.EQ.O) GO TO 260 
GAS1R=~4 . 0001 1 E— 0 1 
GAS2R=2.24484E-01*Z 

GAS3R=(-7.4755E-03+4.86322E-03*Z+6.74265E-04*Y)*Y 

GAS4R=-2.79517E-02*Z*Z 

GAS2E=GAS2/Z 

GAS3E=2.43161E-03*Y*Y 

GAS4E=(-2 . 59512E01— 5 . 59034E-02*Y+3 . 68994E00*Z)*Z 
GAMM R=GAS 1 R+GAS2 R+GAS 3 R+GAS 4R 
GAMME=GAS2E+GAS3E+GAS4E 
GO TO 260 

150 IF(Z.GT.0.65E00) GO TO 160 

GAMM=1 . 398E00 
GO TO 250 

160 IF (Z.GT.1.5E00) GO TO 170 

GAS 1=1 . 39123E00-4 . 08321E-03*Y 

GAS2=( 1 . 42545E-02+1 . 41769E-02*Y)*Z 

GAS3=(2 . 57225E-04+6 . 52912E-04*Z+8 . 46912E~05*Y)*Y*Y 

GAS4=( 6 . 2555E-02-7 . 83637E-03*Y-9 . 78720E-02*Z)*Z*Z 

GAS5=5. 80955-1. 82302E-01*Y 

GAS6=(-9 . 62396E00+1 . 79619E~01*Y)*Z 

GAS7=(-2 . 30518E-02+1 . 18720E-02*Z-3 . 35499E-04*Y)*Y*Y 

GAS8=( 5 . 27047E00-3 . 65507E-02*Y-9 . 19897E-01*Z)*Z*Z 

GAS9=(-10 . 0E00*Z+14.2E00) 

IF(KFLAG.EQ.O) GO TO 240 
GASlR=-4. 08321E-03 
GAS2R=1 . 41769E-02*Z 

GAS3R=( 5 . 1445E-04+1 . 305824E-03*Z+2 . 540736E-04*Y)*Y 
GAS4R=-7 . 83637E-03*Z*Z 
GAS5R=-1 .82302E-01 
GAS6R=1 . 79619E— 01*Z 

GAS7R=( -4 . 61036E-02+2 . 3744E-02*Z-1 . 006497E-03*Y)*Y 

GAS8R=-3.65507E-02*Z*Z 

GAS9R=0.0E00 

GAS2E=GAS2/Z 

GAS3E=6. 52912E-04*Y*Y 

GAS4E=( 1 . 2511E-01-1 . 567274E-02*Y-2 . 936 16E-0 1*Z ) *Z 

GAS6E=GAS6/Z 

GAS7E=1.1872E-02*Y*Y 

GAS8E=( 1 . 054094E01-7 . 31014E-02*Y-2 . 759691E00*Z)*Z 
GAS9E=-10 . 0E00 
GO TO 240 

170 IF (Z.GT.2.22E00) GO TO 180 
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GAS 1=— 1 . 20784E00-2 . 57909E-01*Y 
GAS2=( 5 . 02307E00+2 . 87201E-01*Y)*Z 
GAS3=(-9 . 95577E-03+5 . 23524E-03*Z-1 . 45574E~04*Y)*Y*Y 
GAS4=( -3 . 20619E0O-7 . 50405E-02*Y+6 . 51564E-01*Z)*Z*Z 
GAS 5= -6 . 62841E00+2 . 77112E-02*Y 
GAS 6= (7 . 30762E00-7 . 68230E-02*Y)*Z 
GAS7=(7 . 19421E-03-3 . 62463E-03*Z+1.62777E-04*Y)*Y*Y 
GAS8=("2 . 33161E00+3 . 04767E-02*Y+1 . 66856E-01*Z)*Z*Z 
GAS9=( 1 . 255324E02+2 . 015335E00*Y-6 . 39Q747E01*Z- 
*6. 515225E-01*Y*Z) 

IF(KFLAG.EQ.O) GO TO 240 
GASlR=-2 . 57909E-O1 
GAS2R=2.87201E-01*Z 

GAS3R=( -1 . 991154E-02+1 . 047048E-02*Z-4 . 36722E-04*Y)*Y 

GAS4R=-7.50405E-02*Z*Z 

GAS5R=2 . 77112E-02 

GAS6R=-7.6823E-02*Z 

GAS7R=( 1 . 438842E-02-7 . 24926E-03*Z+4 . 88331E-04*Y)*Y 

GAS8R=3.04767E-02*Z*Z 

GAS9R=2 . 015335E00-6 . 515225E-01*Z 

GAS2E=GAS2/Z 

GAS3E=5 . 23524E-03*Y*Y 

GAS4E=(-6.41238E00-1.50081E-Ol*Y+1.954692E00*Z)*Z 

GAS6E=GAS6/Z 

GAS7E=-3.62463E-03*Y*Y 

GAS8E=( -4 . 66322E00+6 . 09534E-02*Y+5 . 00568E-01*Z) *Z 
GAS9E=-6 . 390747E01-6 . 515225E-01*Y 
GO TO 240 

180 IF (Z.GT.2.95) GO TO 190 

GASl=-2 . 26460EOO-7 . 82263E-02*Y 
GAS2=( 4 . 90497E00+7 . 18096E-02*Y)*Z 
GAS3=( ~3 . 06443E-03+1 . 74209E-03*Z+2 . 84214E-05*Y)*Y*Y 
GAS4=(-2 . 24750E00-1 . 31641E-02*Y+3 . 33658E-01*Z)*Z*Z 
GAS5=-1 . 47904E01-1 . 76627E-01*Y 
GAS6=( 1 . 35036E01+8 . 77280E-02*Y)*Z 
GAS7=(-2 . 13327E-03+7 . 15487E-04*Z+7 . 30928E~05*Y)*Y*Y 
GAS8=(-3.95372EOO-8.96151E-03*Y+3.63229E-01*Z)*Z*Z 
GAS9=( 1 . 788542E02+6 . 317894E00*Y-6 . 756741EQ1*Z- 
*2 . 460060E00*Y*Z) 

IF(KFLAG.EQ.O) GO TO 240 
GASlR=-7 .82263E-02 
GAS2R=7.18096E-02*Z 

GAS3R=( -6 . 12886E-03+3 . 48418E-03*Z+8 . 52642E-05*Y)*Y 

GAS4R=-1.31641E-02*Z*Z 

GAS5R=-1 . 76627E-01 

GAS6R=8. 7728E~02*Z 

GAS7R=(-4 . 26654E-03+1 . 430974E-03*Z+2 . 192784E-04*Y)*Y 

GAS8R=-8.96151E-03*Z*Z 

GAS9R=6 . 317894E00-2 . 46006E00*Z 

GAS2E=GAS2/Z 
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GAS3E=1 . 74209E-03*Y*Y 

GAS4E=(-4 . 495E00-2 . 63282E-02*Y+1 . 000974E00*Z)*Z 

GAS6E=GAS6/Z 

GAS7E=7 . 15487E-04*Y*Y 

GAS8E=(-7 . 90744E00-1. 792302E-02*Y+1 . 089687E00*Z)*Z 
GAS9E=-6 . 756741E01— 2 . 460O6E00*Y 
GO TO 240 

190 IF(Z.GT.3.4E0G) WRITE(6 , 1000) R,E 

GAS1=-1. 66904E01-2 . 58318E-01*Y 
GAS2=( 1 . 78350E01+1 . 54898E-01*Y)*Z 
GAS3=(-9 . 71263E-03+3 . 97740E-03*Z+9 .04300E-05*Y)*Y*Y 
GAS4=(-5 . 94108E00-2 . 01335E-02*Y+6 . 60432E~01*Z)*Z*Z 
GAS 5=8 . 54690E01+1 . 17554E01*Y 
GAS6=(-7 . 21760E01-7 . 15723E00*Y)*Z 
GAS7=(-4. 16150E-02+1 . 38147E-02*Z+5.45184E-04*Y)*Y*Y 
GAS8=( 2 . 01758EQ1+1 . G8990E00*Y-1 . 86438E00*Z)*Z*Z 
GAS9=(2 .883262E02+1 . 248536E01*Y-8.816985E01*Z- 
*3 . 720309E00*Y*Z) 

IF(KFLAG.EQ.O) GO TO 240 
GASlR=-2 . 583 18E— 0 1 
GAS2R=1.54898E-01*Z 

GAS3R=(-1 . 942526E-02+7 . 9548E-03*Z+2 . 7129E-04*Y)*Y 

GAS4R=-2.01335E-02*Z*Z 

GAS5R=1 . 17554E01 

GAS6R=-7 . 15723E00*Z 

GAS7R=(-8 . 323E-02+2 . 76294E-02*Z+1 . 635552E-03*Y)*Y 

GAS8R=1 . 0899E00*Z*Z 

GAS9R=1 . 248536E01-3 . 720309E00*Z 

GAS2E=GAS2/Z 

GAS3E=3 . 9774E-03*Y*Y 

GAS4E=( — 1 . 1882 16E0 1-4. 0267E~02*Y+1 . 981296E00*Z)*Z 

GAS6E=GAS6/Z 

GAS7E=1.38147E-02*Y*Y 

GAS8E=(4 .03516E01+2 . 1798E00*Y-5 . 59314E00*Z)*Z 
GAS9E=-8 . 816985E01-3 . 720309E00*Y 
GO TO 240 

200 IF(Z .GT . 0 . 65E00) GO TO 210 
GAMM=1 . 3988E00 
GO TO 250 

210 IF(Z .GT. 1 . 7E00) GO TO 220 

GAS 1=1 . 37062E00+1 . 29673E-02*Y 

GAS2=( 1 . 1 1418E— 01— 3 . 26912E-02*Y)*Z 

GAS3=( 1 . 06869E-03-2 . 00286E-03*Z+2. 38305E-04*Y)*Y*Y 

GAS4=(-1 . 06133E-01+1 ,90251E-02*Y+3 . 02210E-03*Z)*Z*Z 

GAMM=GAS1+GAS2+GAS3+GAS4 

IF(KFLAG.EQ.O) GO TO 260 

GAS1R=1 . 29673E-02 

GAS2R=-3.26912E-02*Z 

GAS3R=(2.13738E-03-4.00572E-03*Z+7.14915E-04*Y)*Y 

GAS4R=1.90251E-02*Z*Z 


73 



GAS2E=GAS2/Z 

GAS3E=-2.00286E-03*Y*Y 

GAS4E=( -2 . 12266E-01+3 . 80502E-02*Y+9 . 0663E-03*Z)*Z 
GAMMR=GAS 1 R+GAS2 R+GAS3 R+GAS4R 
GAMME=GAS 2 E+GAS 3 E+GAS4E 
GO TO 260 

220 IF(Z.GT.2.35) GO TO 230 

GAS 1=3 . 43846E-02-2 . 33584E-01*Y 

GAS2=( 2 . 85574E00+2 . 59787E-01*Y)*Z 

GAS3=( -10 . 89927E-03+4 . 23659E-03*Z+3 . 85712E-04*Y)*Y*Y 

GAS4=(-1 . 94785E00-6 . 73865E-02*Y+4 . 08518E-01*Z)*Z*Z 

GAS5=-4 . 20569E00+1 . 33139E-01*Y 

GAS6=(4. 51236E00-1 . 66341 E-01*Y)*Z 

GAS7=( 1 . 67787E-03-1 . 10022E-03*Z+3 . 06676E~04*Y)*Y*Y 

GAS8=(-1 . 35516E00+4 . 91716E-02*Y+7 . 52509E-Q2*Z)*Z*Z 

GAS9=( 1 . 757042E02-2 . 163278E00*Y-8 . 833702E01*Z+ 

*1 .897543E00*Y*Z) 

IF(KFLAG.EQ.O) GO TO 240 
GASlR=-2 . 33584E-01 
GAS2R=2.59787E-01*Z 

GAS3R=(-2 . 179854E-02+8 . 47318E-03*Z+1 . 157136E-03*Y)*Y 

GAS4R=-6.73865E-02*Z*Z 

GAS5R=1 . 33139E— 01 

GAS6R=-1 . 66341 E— 0 1 *Z 

GAS7R=( 3 . 35574E-03-2 . 20044E-03*Z+9 . 20028E-04*Y)*Y 

GAS8R=4.91716E-02*Z*Z 

GAS9R=-2 . 163278E00+1 . 897543E00*Z 

GAS2E=GAS2/Z 

RAS3F=4 ?36F9F-03*Y*Y 

GAS4E=(-3 . 8957E00-1 . 34773E-01*Y+1 . 225554E00*Z)*Z 

GAS6E=GAS6/Z 

GAS7E=-1.10022E-03*Y*Y 

GAS8E=( -2 . 71032E00+9 . 83432E-02*Y+2 . 257527E-01*Z)*Z 
GAS9E=-8 . 833702E01+1 . 897543E00*Y 
GO TO 240 

230 IF(Z.GT.2.9E00) WRITE(6,1000) R,E 

GAS1=-1 . 70633E00-1 . 48403E-01*Y 
GAS2=(4. 23104E00+1 . 37290E-01*Y)*Z 
GAS3=(-9.10934E-03+3.85707E-03*Z+2.69026E-04*Y)*Y*Y 
GAS4=( -1 . 97292E00-2 . 81830E-02*Y+2 . 95882E~01*Z)*Z*Z 
GAS5=3 . 41580E01-1 .89972E01*Y 
GAS6=( -4 . 0858E01+1 . 30321E01*Y)*Z 
GAS7=(-8. 01272E-01+2 . 75121E-01*Z-1 . 77969E-04*Y)*Y*Y 
GAS8=( 1 . 60826E01-2 . 23386E00*Y-2 . 08853E00*Z)*Z*Z 
GAS9=( 2 . 561323E02+1 . 737089E02*Y-9 . 058890E01*Z- 
*5 .838803E01*Y*Z) 

IF(GAS9 .GT . 30 . EOO) GAS9=30 . EOO 
I F( GAS9.LT. -30. EOO) GAS9=-30.E00 
IF(KFLAG.EQ.O) GO TO 240 
GAS1R=-1 . 48403E-01 
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GAS2R=1.3729E-01*Z 

GAS3R=( -1 . 821868E-02+7 ,71414E-03*Z+8 . 07078E-04*Y)*Y 

GAS4R=-2.8183E-02*Z*Z 

GAS5R=-1 .89972E01 

GAS6R=1 . 30321 E0 1*Z 

GAS7R=(-1 . 602544E00+5 . 50242E-0 1 *Z-5 . 33907E-04*Y)*Y 

GAS8R=-2.23386E00*Z*Z 

GAS9R=1 . 737089E02-5 . 838803E01*Z 

GAS2E=GAS2/Z 

GAS3E=3.85707E-03*Y*Y 

GAS4E=(-3 . 94584E00-5 . 6366E-02*Y+8 .87646E-01*Z)*Z 

GAS6E=GAS6/Z 

GAS7E=2 . 75121E-Q1*Y*Y 

GAS8E=( 3 . 21652E01-4 . 46772E00*Y-6 . 26559E00*Z)*Z 
GAS9E=-9 . 05889E01-5 . 8388G3E01*Y 
240 GAS9=EXP(GAS9) 

GAMM=GASl+GAS2+GAS3+GAS4+(GAS5+GAS6+GAS7+GAS8)/( 1 . +GAS9) 
IF(KFLAG.EQ.O) GO TO 260 

GAMMR=GAS1R+GAS2R+GAS3R+GAS4R+(GAS5R+GAS6R+GAS7R+GAS8R)/ 

* ( 1 . +GAS9 ) -( GAS5+GAS6+GAS7+GAS8 ) *GAS9R*GAS9/( ( 1 . +GAS9 ) **2 ) 
GAMME=GAS2E+GAS3E+GAS4E+(GAS6E+GAS7E+GAS8E)/( 1 .+GAS9) 
*-(GAS5+GAS6+GAS7+GAS8)*GAS9E*GAS9/((l.+GAS9)**2) 

GO TO 260 
250 CONTINUE 

IF(KFLAG.EQ.O) GO TO 260 
GAMMR=0 . OEOO 
GAMME=0 . OEOO 
260 P=(GAMM-1)*E*R 

IF(KFLAG.EQ.O) GO TO 270 
GAMMR=GAMMR/2 . 302585E00 
GAMME=GAMME/2 . 302585E00 
ASQ=E*( (GAMM-1 . EOO ) * ( GAMM+GAMME ) +GAMMR) 

A=SQRT(ASQ) 

270 I F( I FLAG) 10,280,290 

280 IF(JFLAG) 30,40,10 

290 IF(JFLAG) 60,70,10 

300 X=ALOG10(P/P0) 

Y=ALOG10( R/RO) 

Z1=X-Y 

IF (Y.GT.-0.5E00) GO TO 400 
IF(Y.GT.-4.5E0G) GO TO 350 
IF(Z1 .GT . 0 . 25E00) GO TO 310 
T=P/(GASCON*R) 

RETURN 

310 IF(Z1 . GT . 0 . 95E00) GO TO 320 

GAS 1=1 . 44824E-01+1 . 36744E-Q2*Y 

GAS2=( 1 . 17099E-01-8. 22299E-02*Y)*Z1 

GAS3=(-6 . 75303E— 04— 1 . 47314E-03*Zl-7 . 90851E~05*Y)*Y*Y 

GAS4=( 1 . 3937E00+6 . 83066E-02*Y-6 . 65673E-Q1*Z1)*Z1*Z1 

TN0N=GAS1+GAS2+GAS3+GAS4 
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GAS 5=9 . 07979E00+1 . 01308E00*Y 
GAS6=(-2 . 29428E01-1 . 52122E00*Y)*Z1 
GAS7=(3. 78390E-02-2 . 63115E-02*Zl+5 . 46402E-04*Y)*Y*Y 
GAS8=( 1 . 95657E01+5 . 73839E-01*Y-5 . 63057E00*Z1)*Z1*Z1 
GAS9=EXP(7.619803E01-1.501155E01*Y-6.770845E01*Z1 
*+l . 273147E01*Y*Z1) 

GO TO 440 

380 IF (Z1.GT.2.00E00) GO TO 390 
GAS1=-1 . 1 3598E0 1— 1 . 02049E00*Y 
GAS2=(2 . 22793E01+1 . 24038E00*Y)*Z1 
GAS3=(-3 . 10771E-02+1 . 92551E-02*Z1~2 . 69140E-04*Y)*Y*Y 
GAS4=(-1 . 31512E01-3 . 62875E-01*Y+2 . 64544EOO*Z1)*Z1*Z1 
GAS 5=8 . 72852E00+1 . 27564E00*Y 
GAS6=(-1 . 79172E01— 1 . 52051E00*Y)*Z1 
GAS7=(4 . 91264E-02-2 . 81731 E-02*Zl+5 . 23383E-04*Y)*Y*Y 
GAS8=( 1 . 16719E01+4 . 45413E-Gl*Y-2 . 45584E00*Z1)*Z1*Z1 
GAS9=EXP( 1 . 84792E02+9 . 583443E00*Y-1 . 020835E02*Z1 
*-4. 166727E00*Y*Z1) 

GO TO 440 

390 GAS1=-1. 76079E0 1—1 . 26579E00*Y 

GAS2=(2 . 48544E01+1 . 09442E00*Y)*Z1 
GAS3=(-3 . 65534E-02+1 . 54346E-02*Z1~4 . 59822E-04*Y)*Y*Y 
GAS4=( — 1 . 08166E01-2 . 27803E-01*Y+1 .60641E00*Z1)*Z1*Z1 
GAS 5=2 . 60669E01+2 . 31791E00*Y 
GAS6=(-3 . 22433E01— 1 . 82645E00*Y)*Z1 
GAS7=(4 . 94621E— 02— 1 . 85542E-02*Zl+5 . 04815E-04*Y)*Y*Y 
GAS8=( 1 . 33829E01+3 . 59744E-01*Y~1 .86517E00*Z1 )*Z1*Z1 
GAS9=EXP( 3 . 093755E02+1 . 875018E01*Y-1 . 375004E02*Z1 
*-8.333418E00*Y*Zl) 

GO TO 440 

400 IF (Z1.GT.0.25E00) GO TO 410 
T=P/(GASCON*R) 

RETURN 

410 IF (Z1.GT.0.95E00) GO TO 420 

GASl=-2 . 94081E-03+5 . 73915E-04*Y 

GAS2=(9 .88883E-01-3 . 71241E-03*Y)*Z1 

GAS3=( 1 . 12387E-04-3 . 76528E-04*Z1+1 . 76192E-05*Y)*Y*Y 

GAS4=(2.86656E-02+4.56059E-03*Y-1.99498E-02*Zl)*Zl*Zl 

TN0N=GAS1+GAS2+GAS3+GAS4 

GO TO 450 

420 IF (Z1.GT.1.45E00) GO TO 430 
GAS 1=1 . 32396E00+8 . 52771E-02*Y 
GAS2=(-3 . 24257E00-2 . 00937E-01*Y)*Z1 
GAS3=( 5 . 68146E-03-6 . 85856E-03*Z1+1 . 98366E-04*Y)*Y*Y 
GAS4=(4. 53823E00+1 . 18123E-01*Y~1 .6246E00*Z1)*Z1*Z1 
GAS5=-5 . 26673E— 01— 1 . 58691E~01*Y 
GAS6=(2 . 61600E00+3 . 16356E-01*Y)*Z1 
GAS7=(-1 . 90755E-02+1 . 70124E-02*Zl-5 . 58398E-04*Y)*Y*Y 
GAS8=(-3 . 3793E00— 1 . 52212E-0 1*Y+1 . 30757E00*Z1)*Z1*Z1 
GAS9=EXP( 1 . 442206E02-2 . 544727E01*Y-1 . 277055E02*Z1 
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SUBROUTINE TGAS2(E,R,S) 

INPUTS FOR SUBROUTINE: 
E=INTERNAL ENERGY IN (M/SEC)**2 
R=DENSITY IN KG/M**3 

OUTPUT: 

S=ENTR0PY IN (M/SEC)**2/K 


DATA EO , RO , GASCON/78408 . 4E00 , 1 . 292E00 , 287 . 06E00/ 
RRATIO=R/RO 
ERATI0=E/E0 
Y=ALOG10(RRATIO) 

Z=AL0G10( ERATIO) 

I F(ABS(Y+4 . 5E00) . LT . 2 . 5E-02) GO TO 10 
IF(ABS( Y+0 . 5E00) . LT . 0 . 5E-02) GO TO 40 
IFLAG=-1 
GO TO 80 
10 IFLAG=0 
RSAVE=R 
YM=Y 

Y=-4 . 5E00+2 . 5E-02 
YHIGH=Y 
R=( 10 . **Y)*R0 
JFLAG=-1 
GO TO 80 
20 SHIGH=S 

Y=-4 . 5E00-2 . 5E-02 
YLOW=Y 

R=( 10 . **Y)*R0 
JFLAG=0 
GO TO 80 
30 SLOW=S 

GO TO 70 
40 IFLAG=1 

RSAVE=R 
YM=Y 

Y=-0 . 5E00+0 . 5E-02 
YHIGH=Y 
R=( 10 . **Y)*RG 
JFLAG=-1 
GO TO 80 
50 SHIGH=S 

Y=-0 . 5EOO-0 . 5E-02 
YLOW=Y 

R=( 10. **Y)*RQ 
JFLAG=0 
GO TO 80 
60 SLOW=S 
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S=SLOW+(SHIGH-SLOW)/(YHIGH-YLOW)*(YM-YLOW) 

R=RSAVE 
RETURN 

80 CONTINUE 

I F( Z . LE.O. 65E00) GO TO 110 

IF(Y.GT.-4. 5E00) GO TO 90 

IF(Z.GT.3.69E0G) WRITE(6,1000) R,E 

GASl=-9 . 91081E-01-5 . 00277E00*Y 

GAS2=( 5 . 46521EQ1+5 . 10144E00*Y)*Z 

GAS3=( 1 . 76206E-02+2 . 12002E-02*Z+1 . 76358E-03*Y)*Y*Y 

GAS4=(-2.97001E01-1.84915E00*Y+5.87892E00*Z)*Z*Z 

GO TO 120 

90 IF(Y.GT.-0.50E00) GO TO 100 

IF(Z.GT.3.4E00) WRITE(6, 1000) R,E 

GAS 1=1 . 0836E01-4 . 55524E00*Y 

GAS2=(2 . 96473E01+3 . 90851E00*Y)*Z 

GAS3=(-2 . 05732E-03+3 . 65982E-02*Z+5 . 23821E-03*Y)*Y*Y 

GAS4=(-1 .67001E01-1 .44623E00*Y+3. 98307E00*Z)*Z*Z 

GO TO 120 

100 IF(Z.GT.3.0E0G) WRITE(6 , 1000) R,E 

GAS 1=2 . 01858E01-3 . 13458E00*Y 
GAS2=( 1 . 03619E01+1 .87767E00*Y)*Z 

GAS3=(-1.72922E-01+1.12174E-01*Z+1.28626E-02*Y)*Y*Y 
GAS4=( -5 . 43557E00-8 . 71048E-01*Y+2 . 01789E00*Z)*Z*Z 
GO TO 120 

110 DELTZ=Z-0 . 4E00 

DELTS=(2 . 5EOO*DELTZ-Y)*GASCON*2 . 302585E00 
S=6779 . 2004E00+DELTS 
GO TO 130 

120 SN0N=GAS1+GAS2+GAS3+GAS4 
S=GASCON*SNON 

130 I F( I FLAG) 160,140,150 

140 IF(JFLAG) 20,30,160 

150 IF(JFLAG) 50,60,160 

160 CONTINUE 

1000 FORMAT ( /20X , 48HWARNI NG ! OUTSIDE VALIDITY RANGE OF CURVE FIT 
*,/,20X,5HRB0 =, 1PE15.8,5X,3HE =,1PE15.8,/) 

RETURN 

END 
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SUBROUTINE TGAS3 (P,RH0,T) 

INPUTS FOR SUBROUTINE : 

P= PRESSURE, IN NEWT0NS/M**2. 
RH0= DENSITY, IN KG/M**3. 

OUTPUT : 

T= TEMPERATURE, IN KELVIN. 


DATA RO , PO, TO, GASCON/1 . 292E00 , 1 .0133E05 , 273 . 15E00 , 287 . 06E00/ 
Y=ALOG10( RHO/RO) 

X=AL0G10(P/P0) 

I F(A3S( Y+4 . 5E00) . LT . 2 . 5E-02) GO TO 20 
IF(ABS(Y+0.5E00) .LT.5.0E-03) GO TO 50 
IFLAG=-1 
GO TO 90 
10 RETURN 

20 I FLAG=0 

RSAVE=RHO 

YM=Y 

Y=-4 . 5E00+2 . 5E-02 
YHIGH=Y 

RHO=( 10 . **Y)*R0 
JFLAG=-1 
GO TO 90 
30 THIGH=T 

Y=-4 . 5E00-2 . 5E-02 
YL0W=Y 

RHO=(10.**Y)*R0 
JFLAG=0 
GO TO 90 
40 TLOW=T 

GO TO 80 
50 IFLAG=1 

RSAVE=RH0 

YM=Y 

Y=-0 . 5E00+0 . 5E-02 
YHIGH=Y 

RHO=( 10 . **Y)*R0 
JFLAG=-1 
GO TO 90 
60 THIGH=T 

Y=-0 . 5E00-0 . 5E-02 
YLOW=Y 

RHO=( 10 . **Y)*R0 
JFLAG=0 
GO TO 90 
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*+2 . 026986E00*Y*Z1) 

GO TO 230 

140 IF (Z1.GT.0.25E00) GO TO 150 
T=P/(RHO*GASCON) 

GO TO 250 

150 IF (Z1.GT.0.95E00) GO TO 160 
GAS 1=2 . 03910E-02+7 . 67310E-03*Y 
GAS2=(8. 48581E— 01— 2 . 93086E-02*Y)*Z1 
GAS3=(8. 40269E-04-1 . 47701E-03*Zl+3 . 13687E-05*Y)*Y*Y 
GAS4=(2 . 67251E-01+2. 37262E-02*Y-1 . 41973E-01*Z1)*Z1*Z1 
TN0N=GAS1+GA$2+GAS3+GAS4 
GO TO 240 

160 IF (Z1.GT.1.45E00) GO TO 170 
GASl=-5 . 12404E00-2 . 84740E-01*Y 
GAS2=( 1. 54532E01+4 . 52475E-0 1*Y ) *Z 1 
GAS3=(-1 . 2288 1 E-02+8. 56845E~03*Z 1-3 . 25256E-04*Y)*Y*Y 
GAS4=(-1 . 35181E01-1 . 68725E~01*Y+4. 18451E00*Z1)*Z1*Z1 
GAS 5=7 . 52564E00+8 . 35238E-01*Y 
GAS6=(-1 . 95558E01-1 . 23393E00*Y)*Z1 
GAS7=(3. 34510E-02-2 . 34269E-02*Z1+4.81788E-04*Y)*Y*Y 
GAS8=( 1 . 71779E01+4 . 54628E-01*Y-5 . 09936E00*Z1)*Z1*Z1 
GAS9=EXP(6 . 148442E01-1 .828123E01*Y-5 . 468755E01*Z1 
*+l . 562500E01*Y*Z1) 

GO TO 230 

170 IF (Z1.GT.2.05EO0) GO TO 180 
GAS1=-1 . 23779E0 1— 1 . 14728E00*Y 
GAS2=(2 . 41382E01+1 . 38957E00*Y)*Z1 
GAS3=(-3 . 63693E-02+2 . 24265E-02*Zl-3 . 23888E-04*Y)*Y*Y 
GAS4=(-1 . 42844E01-4 . 06553E-01*Y+2 .87620E00*Z1)*Z1*Z1 
GAS5=4 . 40782E00+1 . 33046E00*Y 
GAS6=(-1 . 15405E01-1 . 59892E00*Y)*Z1 
GAS7=(5 . 30580E-02-3 . 10376E-02*Zl+4. 77650E-04*Y)*Y*Y 
GAS8=(8. 57309E00+4. 71274E-01*Y-1 . 96233E00*Z1)*Z1*Z1 
GAS9=EXP( 1 . 4075E02-6 . 499992E00*Y-7 . 75E01*Zl+5 . 0E00*Y*Z1) 
GO TO 230 

180 IF (Z1.GT.2.50E00) WRITE (6,1000) RHO,P 
GAS1=-1 . 27244E0 1-1 . 66684E00*Y 
GAS2=( 1 . 72708E01+1 . 45307 E00*Y)*Z1 

GAS3=(-3.64515E-02+1.90463E-02*Zl+4.80787E-04*Y)*Y*Y 
GAS4=(-6.97208E00-3.04323E-01*Y+9.67524E-01*Z1)*Z1*Z1 
GAS5=7 . 71330E00+5 . 08340E-01*Y 
GAS6=(-9 .82110E00-4. 49138E-01*Y)*Z1 
GAS7=(-9.41787E-04-2.40293E-03*Zl-8.28450E-04*Y)*Y*Y 
GAS8=(4 . 16530E00+9 . 63923E-02*Y-5 . 88807E-01*Z1)*Z1*Z1 
GAS9=EXP(-1 . 092654E03-3 . 05312E02*Y+4 . 656243E02*Z1+ 

*1 . 312498E02*Y*Z1) 

GO TO 230 

190 IF (Z1.GT.0.25E00) GO TO 200 
T=P/( RHO*GASCON ) 

GO TO 250 
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200 IF (Z1.GT.1.00E00) GO TO 210 

GA$1=-1 . 54141E-03+6 . 58337E-04*Y 

GAS2=(9 .82201 E-0 1-3 .85028E~03*Y)*Z1 

GAS3=(1.23111E-04-4.08210E-04*Zl+2.13592E-05*Y)*Y*Y 

GAS4=( 3 . 77441E-02+4 . 56963E-03*Y-2 . 35172E-02*Z1)*Z1*Z1 

TN0N=GAS1+GAS2+GAS3+GAS4 

GO TO 240 

210 IF (Z1.GT.1.45E00) GO TO 220 

GAS1=8.06492E-01+9.91293E-02*Y 
GAS2=(-1 . 70742E00-2 . 28264E-01*Y)*Z1 
GAS3=( 5 . 03500E-03-6 . 13927E~03*Z1+1 . 69824E-04*Y)*Y*Y 
GAS4=(3 . 02351E00+1 . 31574E-01*Y-1 . 12755E00*Z1)*Z1*Z1 
GAS5=-1 . 17930E-01-2 . 12207E-01*Y 
GAS6=( 1 . 36524E00+4 . 05886E-01*Y)*Z1 
GAS7=(-1 .88260E-02+1 ,65486E-02*Zl-5. 11400E-04*Y)*Y*Y 
GAS8=(-2. 1G926E00-1 .89881 E-01*Y+8. 79806E-01*Z1)*Z1*Z1 
GAS9=EXP( 1 . 959604E02-4 . 269391 EO 1*Y-1 . 73493 1 E02*Z 1+ 
*3.762898E01*Y*Z1) 

GO TO 230 

220 IF (Z1.GT.2.3E00) WRITE (6,1000) RHO,P 
GAS1=-1 . 66249E00-8. 91113E-02*Y 
GAS2=(4. 11648E00+8 . 78093E-02*Y)*Z1 
GAS3=(-3 . 09742E-03+1. 99879E-03*Zl+6 . 85472E-05*Y)*Y*Y 
GAS4=(-1.84445E00-7.50324E-03*Y+3.05784E-01*Z1)*Z1*Z1 
GAS5=1 . 11555E01+1 . 32100E00*Y 
GAS6=(-1 . 7 1236E0 1— 1 . 29190E00*Y)*Z1 
GAS7=(6.28124E-02-3.07949E-02*Zl+1.57743E-03*Y)*Y*Y 
GAS8=(8. 63804E00+3 . 07809E-01*Y-1 . 42634E00*Z1)*Z1*Z1 
GAS9=EXP( 1 . 330611E02+8 . 979635E00*Y-7 . 265298E01*Z1 
*-2 . 449009ED0*Y*Z1) 

230 TN0N=GASl+GAS2+GAS3+GAS4+(GAS5+GAS6+GAS7+GAS8)/( 1 . +GAS9) 

240 T=( 10 . **TNON)*TO 

250 IF (I FLAG) 10,260,270 

260 IF (JFLAG) 30,40,10 

270 IF (JFLAG) 60,70,10 

1000 FORMAT ( /20X , 48HWARN I NG ! OUTSIDE OF VALIDITY RANGE OF CURVE FIT 
*,/,20X,5HRH0 =,1PE15.8,5X,3HP =, 1PE15.8,/) 

END 
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SUBROUTINE TGAS4 (P.RHO.H) 


INPUTS FOR SUBROUTINE : 

P = PRESSURE, IN NEWTONS/M. **2. 

RHO = DENSITY, IN KG/M**3. 

OUTPUT : 

H = SPECIFIC ENTHALPY, IN (M/$EC)**2. 


DATA RO , PO/1 . 292E00 , 1 . 0133E05/ 
Y=ALOG10( RHO/RO) 

X=ALOG10(P/P0) 

I F(ABS( Y+4 . 5EOO) . LT . 2 . 5E-02) GO TO 20 
I F(ABS( Y+0 . 5E00) . LT . 5 . 0E-03) GO TO 50 
IFLAG=-1 
GO TO 90 
10 RETURN 

20 IFLAG=0 

RSAVE=RHO 

YM=Y 

Y=-4 . 5E00+2 . 5E-02 
YHIGH=Y 

RHO=( 10 . **Y)*R0 
JFLAG=-1 
GO TO 90 
30 HHIGH=H 

Y=-4 . 5E00-2 . 5E-02 
YLOW=Y 

RHO=( 10 . **Y)*R0 
JFLAG=0 
GO TO 90 
40 HL0W=H 

GO TO 80 
50 IFLAG=1 

RSAVE=RHO 

YM=Y 

Y=-0 . 5E00+0 . 5E-02 
YHIGH=Y 

RHO=( 10 . **Y)*R0 
JFLAG=-1 
GO TO 90 
60 HHIGH=H 

Y=-0 . 5E00-0 . 5E-02 
YLOW=Y 

RHO=( 10. **Y)*R0 
JFLAG=0 
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230 

240 


GAMM=GASl+GAS2+GAS3+GAS4+(GAS5+GA$6+GAS7+GA$8)/( 1.+GAS9) 
H=GAMM/(GAMM-1.0E00)*P/RHO 
IF (I FLAG) 10,250,260 
250 IF (JFLAG) 30,40,10 

260 IF (JFLAG) 60,70,10 

1000 F0RMAT(/20X,48HWARNING! OUTSIDE OF VALIDITY RANGE OF CURVE FIT 
*,/,20X,5HRHO =,1PE15.8,5X,3HP =, 1PE15 .8,/) 

END 
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oooooooooo 


SUBROUTINE TGAS5(P,S,RH0) 


INPUTS FOR SUBROUTINE: 

P=PRESSURE IN NEWT0NS/M**2 . 

S=SPECI FI C ENTROPY IN (M/SEC)**2 ./KELVIN 

OUTPUT: 

RHO=DENSITY IN KG/M**3. 


DATA P0,R0, SO, GASCON/1 .0133E05 , 1 . 292E00 , 6779 . 2E00 , 287 . 06E00/ 
X=ALOG10(P/P0) 

SNON=ALOG10( S/GASCON) 

Z=X-SN0N 

IF (SNON.GE. 1.23E00) GO TO 10 
DELTS=S-S0 

RRATIO=ALOG( P/P0)/1 . 4-DELTS/( 3 . 5*GASC0N) 

RHO=RO*EXP( RRATIO) 

RETURN 

10 IF (SNON.GE. 1.42E00) GO TO 20 

GAS1=-1 . 72119E01+5 . 49354E01*SNON 
GAS2=(-1 . 99776E00+3 . 17884EOO*SNON)*Z 

GAS3=(-4.69831E01-8.66580E-01*Z+1.21069E01*SNON)*SNON*SNON 
GAS4=( 1 . 58567E-01-1 . 03055E-01*SN0N-1 . 52322E-03*Z)*Z*Z 
Y=GAS1+GAS2+GAS3+GAS4 
GO TO 120 

20 IF (SNON.GE. 1.592E00) GO TO 30 

GASl=-2 . 78074E02+6 . 11791E02*SNON 
GAS2=( 1 . 37528E0 1—1 . 92394E01*SNON)*Z 

GAS3=(-4. 42909E02+7 . 10425E00*Z+1 . 05869E02*SNON)*SNON*SNON 
GAS4=(-1 . 97269E-01+1 . 49708E-01*SNON+1 . 19153E-03*Z)*Z*Z 
GAS5=2 . 80393E02-5 . 95834E02*SN0N 
GAS6=(~1 . 75934E01+2.49706E01*SNON)*Z 

GAS7=(4.21767E02-8.85461EOO*Z-9.94515E01*SNON)*SNON*SNON 
GAS8=( 3 .61896E-01-2 . 55458E-01*SN0N-2 . 96892E-03*Z)*Z*Z 
GAS9=EXP( -15 . 0*( X+54 . 179-86 . 947*SNON+33 . 583*SNON*SNON)) 

GO TO 110 

30 IF (SNON.GE. 1.70E00) GO TO 50 
ZM=7 . 5269E00*SN0N-14 . 9366 
IF (Z.GE.ZM) GO TO 40 
GAS 1=1 . 10732E02— 1 . 33968E02*SNON 
GAS2=(3 . 7583E-01+2 . 77887E-01*SNON)*Z 
GAS3=4. 03018E01*SNON*SNON 
GAS4=( 1 . 18506E-01-6 . 98812E-02*SNON)*Z*Z 
GAS 5=-6 .85292 EO 1+8. 23834E01*SNON 
GAS6=(-1 .24942E00+8. 31615E-01*SN0N)*Z 
GAS7=-2 ,4524E01*SNON*SNON 
GAS8=(-1 . 13019E-01+7 . 46719E-02*SNON)*Z*Z 
GAS9=EXP( -2 . 0*( Z+33 . 976*SNON-49 .659)) 
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40 


GO TO 110 

GAS 1=3 . 97149E01~4.80988E01*SNON 
GAS2=( — 1 . 17821E00+1 . 27522EOO*SNON)*Z 
GAS3=1 . 45552E01*SNON*SNON 
GAS4=(-2 . 50279E-01+1 . 58399E-01*SN0N)*Z*Z 
GAS5=-5 . 01859E01+6 . 31564E01*SNON 
GAS6=(2 . 39925E00-1 .47883EOO*SNON)*Z 
GAS7=-1 . 98852E01*SNON*SNON 
GAS8=(4. 66791E-01-3 . 06926E-01*SNON)*Z*Z 
GAS9=EXP(-2.0*(Z+11.16*SNON-18.203)) 

GO TO 110 

50 IF (SNON.GE. 1 .80E00) GO TO 70 
ZM=7 . 5269E00*SN0N-14 . 9366 
IF (Z.GE.ZM) GO TO 60 
GAS1=6. 22639E01-7 .8884E01*SNON 
GAS2=( -2 . 19123E00+1 . 84533EOO*SNON)*Z 
GAS 3=2 . 48413E01*SNON*SNON 
GAS4=(-8 . 76105E-02+5 . 55047E-02*SNON)*Z*Z 
GAS5=-7 . 28463E01+9 . 42375E01*SNON 
GAS6=( 4 . 04979E00-2 . 45467EGO*SNON)*Z 
GAS7=-3 . 02865E01*SNON*SNON 
GAS8=( 1 . 2872E-G1-8. 63721E-02*SNON)*Z*Z 
GAS9=EXP(-2 . 0*(Z+13 . 4576*SN0N-I8. 36)) 

GO TO 110 

60 GASl=-2 . 99578E01+3 . 89998EQ1*SN0N 

GAS2=(2 . 9O256E00-1 . 25088E00*SNON)*Z 

GAS3=-1 . 26641E01*SNON*SNON 

GAS4=(6 . 13867E-02-5 . 17059E-02*SNON)*Z*Z 

GAS5=6 .05677 EO 1-7 . 08307E01*SNON 

GAS6=(-2 . 95622E00+1 . 70504E00*SNON)*Z 

GAS7=2 . 07549E01*SNON*SNON 

GAS8=(-2 . 70962E-01+1 . 80057E-01*SNON)*Z*Z 

GAS9=EXP(-2 . 0*(Z+16 .822*SN0N-29 .139)) 

GO TO 110 

70 IF (SNON.GE. 1.90E00) GO TO 80 

GAS1=6 . 23124E02-6 . 77571E02*SNON 

GAS2=(9 . 1281 1E0 1—9 . 55480E01*SNON)*Z 

GAS3=( 1 . 84603E02+2 . 54274E01*Z)*SNON*SNON 

GAS4=(6 . 97635E-01-3 . 27916E-01*SN0N+4 . 90838E-03*Z)*Z*Z 

GAS5=-6 . 39514E02+6 . 97458E02*SNON 

GAS6=(-1 .00154E02+1 . 06701E02*SN0N)*Z 

GAS7=( — 1 . 90745E02-2 .86323E01*Z)*SNON*SNON 

GAS8=(-4 .81471 E-0 1+2 . 00348E-01*SNON-6 . 43371E~03*Z)*Z*Z 

GAS9=EXP(-2.0*(Z+35.275*SNON-58.624)) 

GO TO 110 

80 IF (SNON.GE. 2. OOEOO) GO TO 90 

GAS 1=1 . 40088E01-1 . 58855E01*SNON 
GAS2=( 1 . 71245E-01+4. 32352E-01*SN0N)*Z 
GAS3=4 . 42836EOO*SNON*SNON 
GAS4=(-9 . 54417E-03+8. 92335E-03*SN0N)*Z*Z 
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GAS5=-1 . 13217E02+1 . 24304E02*SNON 
GAS6=(4 . 41505E00-2 . 48530E00*SN0N)*Z 
GAS7=-3 . 4237OE01*SNON*SNON 
GAS8=(-7 . 84297E-02+1 .21459E-02*SN0N)*Z*Z 
GAS9=EXP( -2 . 00*(Z+20 . 884*SN0N-36 . 54) ) 

GO TO 110 

90 IF (SNON.GT.2. 10) GO TO 100 

GASl=-3 . 62767E01+3 . 8634E01*SN0N 
GAS2=( 1 . 48507E00-3 . 41824E~01*SNON)*Z 
GAS3=-1. 04718E01*SNON*SNON 
GAS4=(-3 . 25437E-02+1 . 22032E-02*SNON)*Z*Z 
GAS5=(6 . 52920E01-6 . 28154E01*SNON) 

GAS6=( 2 . 10906E00-1 . 11759E00*SN0N)*Z 
GAS7=1 . 50026E01*SNON*SNON 
GAS8=(2 . 71716E-01-1. 4431E-01*SNON)*Z*Z 
GAS9=EXP( -5 . 00*(Z+28 . 284*SN0N-53 .185)) 

GO TO 110 

100 GAS 1=1 . 86938E02-1 . 45261E02*SNON 

GAS2=(2 . 23883 EO 1 — 1 . 04974E01*SN0N)*Z 
GAS 3=2 . 6Q284E01*SN0N*SN0N 
GAS4=( 1 .85949E00-8. 98218E-01*SNON)*Z*Z 
GAS5=-1 . 27867E02+8 . 58453E01*SN0N 
GAS6=(-2 . 28808E01+1 . 13034E01*SN0N)*Z 
GAS7=-1 . 1234E01*SNON*SNON 
GAS8=(-1. 56005E00+7 . 79397E-01*SN0N)*Z*Z 
GAS9=EXP( -2 . 00*(Z+15 . 048*SN0N~26 .307)) 

110 Y=GASl+GAS2+GAS3+GAS4+(GAS5+GAS6+GAS7+GAS8)/( 1 . +GAS9) 

120 RHO=( 10 . **Y)*RO 
RETURN 
END 
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ooooooooo 


SUBROUTINE TGAS6(P,S,E) 


INPUTS FOR SUBROUTINE: 

P=PRESSURE IN NEWT0N$/M**2. 

S= ENT ROPY IN (M/SEC)**2. /KELVIN 

OUTPUT: 

E=INTERNAL ENERGY IN (M/SEC)**2. 


DATA PO , EO , SO , GASCON/1 .0133E05 , 78408 . 4E00 , 6779 . 2E00 , 287 . 06E00/ 
X=ALOG10(P/P0) 

SNON=ALOG10( S/GASCON) 

Z=X-SNON 

IF (SN0N.GE.1.23E00) GO TO 10 
DELTS=S-S0 

ERATIO=(ALOG(P/PO)+DELTS/GASCON)/3 . 5E00 
E=2. 5E00*E0*EXP(ERATI0) 

RETURN 

10 IF (SNON.GT.1.4EOO) GO TO 20 

GAS1=4.89511E01-1 . 15989E02*SNON 
GAS2=(-4 . 43026E-01+1 . 49819EOO*SNON)*Z 

GAS3=(8. 70709 EO 1-7 .024 17E-01*Z-2 . 04916E01*SNON)*SNON*SNON 
GAS4=(-5.09747E-02+2.94929E-02*SN0N+5.24685E-04*Z)*Z*Z 
Y=GAS 1+GAS2+GAS3+GAS4 
GO TO 120 

20 IF (SNON.GE. 1 . 592E00) GO TO 30 
GAS 1=1 . 30563E02-2 . 87254E02*SN0N 
GAS2=(-7 . 71873E00+1 . 17605E01*SNON)*Z 

GAS3=(2 .06757E02-4. 31838EOO*Z-4.83472E01*SNON)*SNON*SNON 
GAS4=( 1 . 64299E-01-1 . 21768E-01*SN0N-1 . 26118E-03*Z)*Z*Z 
GAS5=-8 . 4598E01+1 . 74701E02*SN0N 
GAS6=(9 . 52234E00-1 . 32819E01*SN0N)*Z 

GAS7=( — 1 . 20234E02+4 . 6271E00*Z+2 . 75818E01*SNON)*SNON*SNON 
GAS8=(-2 . 27725E-01+1 . 60365E-01*SN0N+1 . 4933E~03*Z)*Z*Z 
GAS9=EXP( -15 . 0*(X+54 . 179-86 . 947*SN0N+33 . 583*SNON*SNON)) 

GO TO 110 

30 IF (SNON.GE. 1.70E00) GO TO 50 
ZM=-1 . 917E00*SNON+0 . 092E00 
IF (Z.GT.ZM) GO TO 40 
GAS1=-1 . 34875E02+1 . 92654E02*SNON 
GAS2=( 1 . 04478E01-6 . 60355E00*SNON)*Z 
GAS3=-6 . 72013E01*SNON*SNON 
GAS4=(4. 14546E-01-2 . 75698E-01*SN0N)*Z*Z 
GAS 5=1 . 077 68E02— 1 . 61662E02*SN0N 
GAS6=(-1 . 01257E01+6 . 4962EOO*SNON)*Z 
GAS7=5 . 93194E01*SNON*SNON 
GAS8=(-4 . 63204E-01+3 . 04448E-01*SN0N)*Z*Z 
GAS9=EXP( -2 . 0*(Z+50 . 04*SNON-74 .698)) 

GO TO 110 
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40 GASl=-3 . 39933E01+4. 02792E01*SNON 

GAS2=(7 . 62742E-01-3 . 79044E-01*SNON)*Z 

GAS3=-1 . 09895E01*SNON*SNON 

GAS4=( -7 . 13412E-02+4 . 22359E-02*SN0N)*Z*Z 

GAS 5=3 . 32182E01-4 . 03847E01*SNON 

GAS6=( 1 . 87601E-01-1 . 07599E-01*SNON)*Z 

GAS 7=1 . 22657E01*SNON*SNON 

GAS8=(-1 . 65537E-01+1 . 02560E-01*SN0N)*Z*Z 

GAS9=EXP( -2 . 0*( Z+14 . 062*SNON-24 . 046) ) 

GO TO 110 

50 IF ( SNON . GE . 1 .80E00) GO TO 70 
ZM=-1.917E00*SNON+0 . 092E00 
IF (Z.GT.ZM) GO TO 60 
GAS1=-1 . 66572E01+1 . 51243E01*SNON 
GAS2=(-4 . 72897E-01+2 . 95848E-01*SN0N)*Z 
GAS3=-2 . 31418EOO*SNON*SNON 
GAS4=(-3 . 61191 E— 02+1 . 7288E-02*SNON)*Z*Z 
GAS5=-3 . 27168E00+9 . 45383EOO*SNON 
GAS6=(2 . 04E00-1 . 1342EOO*SNON)*Z 
GAS7=-4 . 28765EOO*SNON*SNON 
GAS8=( 1 . 5021E-01-8. 27388E-02*SN0N)*Z*Z 
GAS9=EXP(-2 . 0*(Z+16 . 995*SN0N-23 . 354) ) 

GO TO 110 

60 GAS1=-1 . 94882E01+2. 31472E01*SNON 

GAS2=( 1 . 00783E00-5 . 18103E-01*SN0N)*Z 

GAS3=-5 . 92829EOO*SNON*SNON 

GAS4=(9 . 92242E-02-5 . 67928E-02*SN0N)*Z*Z 

GAS5=-2 . 09813E00+1 . 22781EOO*SNON 

GAS6=(-6 . 27531E-01+3 . 71651E-01*SNON)*Z 

GAS7=9 . 84125E-04*SNON*SNON 

GAS8=(-1 . 24024E-01+7 . 65034E-02*SNON)*Z*Z 

GAS9=EXP(-2 . 0*(Z+18. 152*SNON-31 . 096)) 

GO TO 110 

70 IF (SNON.GE. 1.90E00) GO TO 80 

GASl=-6 . 37715E01+5 . 46005E01*SNON 
GAS2=(-8 . 9O073EOO+5 . 02063E00*SNON)*Z 
GAS3=-9 . 67094EOO*SNON*SNON 
GAS4=(-7 . 13493E-01+3 . 97962E-01*SN0N)*Z*Z 
GAS5=5 . 05131E01-3 . 97363E01*SNON 
GAS6=(9. 27248E00-5 . 15576EOO*SNON)*Z 
GAS7=6 . 45632EOO*SNON*SNON 
GAS8=(7.2563E-01-4.03054E-01*SNON)*Z*Z 
GAS9=EXP(-2 . 0*(Z+35 . 275*SN0N-58 . 624)) 

GO TO 110 

80 IF (SNON.GE. 2. OOEOO) GO TO 90 

GAS1=-1 . 75434E01+2 .00307E01*SNON 
GAS2=(6 . 20484E-01-2 . 76756E-01*SNON)*Z 
GAS3=-4 . 76084EOO*SNON*SNON 
GAS4=( 1 . 14921E-02-5 ,80043E-03*SNON)*Z*Z 
GAS5=1 . 10485E01- 1 . 261E01*SNON 
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GA$6=(-6 . 59043E-01+3 . 66266E-01*SNON)*Z 
GAS7=3 . 59417EOO*SNON*SNON 
GAS8=(-6 . 04061E-03+6 . 78382E-03*SNON)*Z*Z 
GAS9=EXP( -2 . 0*(Z+20 . 884*SN0N-36 . 54) ) 

GO TO 110 

90 IF (SNON.GE.2. 10E00) GO TO 100 

GAS1=-1 . 18487E01+1 . 26558E01*SN0N 

GAS2=(6 . 09885E-04+5 . 90396E-02*SN0N)*Z 

GAS3=-2 . 45275EOO*SNON*SNON 

GAS4=(-6 . 70888E-03+5 . 21033E-03*SNON)*Z*Z 

GAS 5= 3 . 39157E00— 3 . 835EOO*SNON 

GAS6=(-2 . 14031 E— 0 1+9 . 90508E-02*SNON)*Z 

GAS7=1 . 070Q8E00*SNON*SNON 

GAS8=(-7 . 87516E-02+3 . 63642E-02*SNON)*Z*Z 

GAS9=EXP( -5 . 0*(Z+28 . 284*SN0N-53 .185)) 

GO TO 110 

100 GASl=-2 . 40656E01+3 . 35589E00*SN0N 

GAS2=(-1 . 14521E01+5 . 62572EOO*SNON)*Z 

GAS3=4 . 93842EOO*SNON*SNON 

GAS4=(-7 . 55608E-01+3 . 69774E-01*SNON)*Z*Z 

GAS5=2 . 23905E01-1 . 90842E00*SN0N 

GAS6=( 1 . 05618E01-5 . 17247EOO*SNON)*Z 

GAS7=-4. 3714EOO*SNON*SNON 

GAS8=(6 . 27818E-01-3 . 12015E-01*SNON)*Z*Z 

GAS9=EXP( -2 . 0*(Z+15 . 048*SNON-26 .307)) 

110 Y=GASl+GAS2+GAS3+GAS4+(GAS5+GAS6+GAS7+GAS8)/( 1 . +GAS9) 

120 E=(10.**Y)*E0 

RETURN 
END 
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SUBROUTINE TGAS7(P,S,A) 

INPUTS FOR SUBROUTINE: 
P=PRESSURE IN NEWT0NS/M**2. 
S=ENTR0PY IN (M/SEC)**2. /KELVIN 

OUTPUT: 

A=SPEED OF SOUND IN M/SEC. 


DATA PO, SO, RO, AO, GASCON/1. 0133E05, 6779. 2E00.1.292E00, 331. 3613E00, 
*287 . Q6E00/ 

X=AL0G10(P/P0) 

SNON=ALOG10( S/GASCON) 

Z=X-SNON 

IF (SNON.GE. 1 . 23E00) GO TO 10 

DELTS=S-S0 

ASQLOG=ALOG( 1 . 4EOO*PO/RO)+(ALOG( P/PO)+DELTS/GASCON)/3 . 5E00 
A=EXP(ASQLOG/2 . 0E00) 

RETURN 

10 IF (SNON.GT. 1.4E00) GO TO 20 

GAS1=-1.38377E-Ol-8.84138EOO*SNON 
GAS2=( 2 . 61050E00-3 . 16535EOO*SNON)*Z 

GAS3=( 1 . 10866E01+9 . 88389E-01*Z-3 . 25761EOO*SNON)*SNON*SNON 
GAS4=( -1 . 00224E-01+6 . 62193E-02*SN0N+8. 20610E-04*Z)*Z*Z 
Y=GAS1+GAS2+GAS3+GAS4 
GO TO 120 

20 IF (SNON.GE. 1.595E00) GO TO 30 
GAS 1=1 . 3 1057E02— 2 . 88847E02*SNON 
GAS2=(-5 . 04887E00+7 . 73862EOO*SNON)*Z 
GAS3=(2 . 10147E02— 2 . 88963E00*Z-5 . 0396E01*SNON)*SNON*SNON 
GAS4=(5.48031E-02-4.39459E-02*SNON-2.10202E-04*Z)*Z*Z 
GAS5=-1 . 33465E02+2 . 84739E02*SNON 
GAS6=( 7 . 57389E00-1 . 07749E01*SNON)*Z 

GAS7=(-2 . 02362E02+3 . 8313E00*Z+4 . 79075ED1*SN0N)*SN0N*SN0N 
GAS8=(-1 . 53453E-01+1 . 08531E-01*SNON+9 . 7931E-04*Z)*Z*Z 
GAS9=EXP( -15 . 0*( X+54 . 179-86 . 947*SNON+33 . 583*SNON*SNON) ) 

GO TO 110 

30 IF (SNON.GE. 1.693E00) GO TO 50 
ZM=-9 . 842E00*SN0N+14 . 19E00 
IF (Z.GT.ZM) GO TO 40 
GASl=-6 . 13548E01+7 . 80742E01*SNON 
GAS2=(2 . 08524E00-1 . 21609E00*SNON)*Z 
GAS3=-2 . 43686E01*SNON*SNON 
GAS4=(8. 77563E-02-5 . 46311E-02*SNON)*Z*Z 
GAS5=2 . 07952E01-2 . 71591E01*SNON 
GAS6=(-7 . 43673E-01+4. 08312E~01*SNON)*Z 
GAS7=8 . 68124EOO*SNON*SNON 
GAS8=( -9 . 3592 E— 02+5 . 32328E-02*SN0N)*Z*Z 
GAS9=EXP( -2 . 0*( Z+38 . 785EOO*SNON-57 . 157E00)) 
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GO TO 110 

40 GAS1=3 . 37056E00-4 . 87016E00*SNON 

GAS2=(-3.85754E-01+2.87192E-01*SN0N)*Z 
GAS3=2 . 02041E00*SN0N*SNON 
GAS4=(-4.63144E-03+8.30832E-03*SN0N)*Z*Z 
Y=GAS1+GAS2+GAS3+GAS4 
GO TO 120 

50 IF (SNON.GE. 1.80) GO TO 70 

ZM=-1 . 917*SN0N+0 . 092 
IF (Z.GT.ZM) GO TO 60 
GASl=-8 . 04927E01+7 . 63739E01*SN0N 
GAS2=(-9 . 381E00+5 . 72104EOO*SNON)*Z 
GAS3=-1 . 63435E01*SNON*SNON 
GAS4=( -7 . 48578E-01+4 . 50043E-01*SNON)*Z*Z 
GAS5=8 . 34054E01-8 . 58837E01*SNON 
GAS6=(4. 84197E00-3 . 11188E0Q*SN0N)*Z 
GAS7=2 . 1 1196E01*SNON*SNON 
GAS8=(2 . 33945E-01-1 . 59099E-01*SNON)*Z*Z 
GAS9=EXP( -2 . 0*( Z+7 . 874*SN0N-7 .569)) 

GO TO 110 

60 GASl=-4 . 73308E00+4 . 69363EOO*SNON 

GAS2=(9 . 43798E-02+3 . 54953E-03*SNON)*Z 
GAS3=-7 . 98293E-01*SNON*SNON 
GAS4=(2 .02561E-02-8. 73036E-03*SNON)*Z*Z 
Y=GAS 1+GAS2+GAS3+GAS4 
GO TO 120 

70 IF (SNON.GE. 1 .900E00) GO TO 80 

GASl=-6 . 60574E02+7 . 38042E02*SNON 

GAS2=( -8 . 77589E01+9 . 7894EQ1*SN0N)*Z 

GAS3=(-2 . 06156E02— 2 . 73753E01*Z)*SNON*SNON 

GAS4=-2.14028E-02*Z*Z 

GAS 5=6 . 65014E02-7 . 43416E02*SN0N 

GAS6=(8.87679E01-9.90508E01*SNON)*Z 

GAS7=(2 .08117E02+2 . 77374E01*Z)*SN0N*SN0N 

GAS8=2 . 81 148E~02*Z*Z 

GAS9=EXP(-2 . 0*(Z+35 . 275*SN0N-58 .624)) 

GO TO 110 

80 IF (SNON.GE. 2. 00E00) GO TO 90 

GASl=-5 . 93554E00-7 . 79929EOO*SNON 

GAS2=(-7 . 23618E00+3 . 31162EOO*SNON)*Z 

GAS3=5 . 06381EOO*SNON*SNON 

GAS4=(-5 . 377735E-01+2 . 46865E-01*SN0N)*Z*Z 

GAS5=3 . 95260E01-2 . 90994E0PSNON 

GAS6=(6 . 24136E00-2 . 70007EOO*SNON)*Z 

GAS7=5 . 42786EOO*SNON*SNON 

GAS8=(5 . 95235E— 0 1—2 . 67771E-01*SN0N)*Z*Z 

GAS9=EXP(-2 , 0*(Z+13 . 6751*SNON-20 .1676)) 

GO TO 110 

90 IF ( SNON.GE. 2. 10E00) GO TO 100 
GAS 1=7 . 02453E01-7 . 34732E0PSN0N 
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GAS2=-5.86844E-01*Z 

GAS3=( 1 . 95548E01+1 . 81521E-01*Z)*SNON*SNON 

GAS4=7.36786E-03*Z*Z 

GAS5=-1 . 69610E01+1 .80169E01*SN0N 

GAS6=-3.65779E-01*Z 

GAS7=(-4.75063E00+1.07583E-01*Z)*SNON*SNON 

GAS8=2 . Q1803E-02*Z*Z 

GAS9=EXP( -5 . 0*( Z+28 . 284*SN0N-53 .185)) 

GO TO 110 

100 GAS1=-1 . 12793E02+1 . 35296E02*SN0N 

GAS2=2 . 74295E00*Z 

GAS3=( —4 . 08171E01— 1 . 29257EOO*Z)*SNON*SNON 

GAS4=-2.25652E-01*Z*Z 

GAS 5=2 . 57057E02-2 . 83199E02*SNON 

GAS6=-5.84656E00*Z 

GAS7=(7 .91574E01+2 . 17312EOO*Z)*SNON*SNON 

GAS8=3.35658E-01*Z*Z 

GAS9=EXP(-2 . 0*(Z+15. 048*SNON-26 .307)) 

110 Y=GASl+GAS2+GAS3+GAS4+(GAS5+GAS6+GAS7+GAS8)/( 1 .+GAS9) 
120 A=( 10 . **Y)*AO 

RETURN 
END 
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